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THE  ATCHAFALAYA  RIVER  DELTA 


ANALYTICAL  ANALYSIS  OF  THE  DEVELOPMENT  OF 
THE  ATCHAFALAYA  RIVER  DELTA 


PART  I;  INTRODUCTION 


Objectives 


1.  The  objective  of  this  research  is  to  quantify  the  dynamic 
interaction  of  the  resources  of  the  Atchafalaya  River,  namely,  the 
sediments  forming  the  delta  and  the  riverflow  that  carries  the 
sediment.  Our  study  is  focused  on  fresh  water  discharging  into  a 
quiescent  bay  and  its  dynamic  response  to  the  forcing  function.  Our 
results  are  developed  for  short-term  predictions  of  the  delta  growth 
in  early  stages. 

2.  The  specific  objectives  of  this  research  are: 

a.  To  apply  the  theory  of  turbulent  jets  in  predicting  the 
short-term  process  of  delta  growth  at  the  river  mouth. 

b.  To  formulate  the  problem  of  river  outlet  freshwater 
discharge  into  a  quiescent  bay  as  a  two-dimensional 
plane  jet. 

c.  To  develop  an  analytical  approach  that  quantifies  the 
areal  and  mass  extent  of  delta  growth  as  influenced  by 
the  river  discharge. 

d.  To  test  the  adequacy  of  the  analytical  technique  based 
on  presently  mapped  bathymetry  and  to  verify  the  result 
with  the  independently  measured  bathmetry  changes  in  the 
bay . 

f.  To  perform  a  sensitivity  analysis  on  various  hydro- 
dynamic  parameters,  and  to  assess  the  relative 
importance  of  river  stage  and  discharge,  channel  con¬ 
figuration,  and  bottom  resistance  as  they  relate  to  the 
river  outlet  sediment  deposition. 

Background  of  Deltaic  Processes 

3.  Coleman  and  Wright  (1975)  discussed  the  various  aspects  of 
interacting  coastal  processes  and  their  effects  on  delta  formation. 
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The  most  important  factors  are  climate,  water  and  sediment  discharge, 
vegetation  and  soil,  geometry  of  river  mouth,  winds  and  nearshore  cur¬ 
rents,  wave  power  and  tidal  regime,  and  bathymetry  of  the  receiving  bay. 
Figure  1  shows  the  major  components  of  a  river  basin,  bay,  and  delta 
system  (Coleman  and  Wright  1975). 

4.  The  river  mouth  is  the  point  at  which  the  fresh  water  leaves 
the  confined  channel  and  spreads  and  mixes  with  ambient  bay  water, 
causing  a  decrease  in  flow  velocity  and  total  momentum,  and  conse¬ 
quently,  the  deposition  of  sediment.  The  river  discharge  depends  on 
the  climatic  and  hydrologic  regime  within  the  drainage  basin.  The 
pattern  of  delta  growth  depends  upon  the  rate  of  sediment  supply  by 
the  river  discharge  and  reworking  of  sediment  by  wave  and  current 
forces  in  the  receiving  bay. 

5 .  The  pattern  of  sediment  deposition  depends  upon  the  proper¬ 
ties  of  the  sediment  and  the  relative  roles  of  three  primary  forces 
(Coleman  1976): 

a.  The  inertial  force  of  river  effluent  and  associated 
turbulent  diffusion. 

b.  The  frictional  force  between  the  river  effluent 
and  the  bed  immediately  seaward  of  the  mouth. 

c.  The  buoyant  force  resulting  from  density  differences 
between  river  effluent  and  ambient  fluids. 

6.  Extensive  observations  and  representative  data  collected  at 
the  mouth  of  the  Mississippi  River  by  Wright  and  Coleman  (1974)  have 
indicated  that  the  relative  roles  of  these  forces  vary  in  space  and 
time,  causing  corresponding  changes  in  the  modes  and  patterns  of 
sediment  transport  and  deposition.  Figure  2  illustrates  the  river 
mouth  mechanisms  and  the  resulting  effluent  plume  and  subaqueous  bar 
(Coleman  1976).  Conditions  for  the  four  cases  illustrated  and  their 
results  are  as  follows: 

a.  Inertial  forces  dominant.  When  riverflow  velocities 
are  high,  depths  immediately  seaward  of  the  mouth  are 
large,  density  differences  are  negligible,  inertial 
forces  are  dominant,  and  the  river  effluent  spreads  and 
diffuses  as  a  turbulent  jet.  Narrow  and  linear  sand¬ 
bars  are  formed  (Figure  2a) . 
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Figure  2.  River  mouth  mechanisms  (Coleman  1976) 


b.  Frictional  forces  dominant.  When  riverflow  velocities 
are  high,  but  the  depths  seaward  of  the  mouth  are 
shallow,  turbulent  flux  penetrates  to  the  entire  water 
column  and  bottom  friction  plays  a  major  role  in  causing 
the  river  effluent  to  be  decelerated  and  expanded  as  a 
fan  shape  plume.  Bifurcating  sandbars  are  established 
(Figure  2b) . 

c.  Buoyant  forces  dominant.  When  the  density  of  ambient 
bay  water  is  much  higher  than  the  density  of  river 
effluent,  then  strong  vertical  density  gradient  exists 
at  the  river  mouth  and  buoyancy  becomes  of  paramount 
importance  in  spreading  both  the  river  effluent  and  sand 
bars  radially  away  from  the  mouth  (Figure  2c). 

d.  Interaction  of  forces.  Various  combinations  of  these 
three  forces  exist  in  modern  deltas.  Interactions 
between  buoyant  and  inertial  forces  are  common  in  many 
modern  rivers  (Figure  2d) . 


Approach  in  Analyzing  the  Atchafalaya  River  Delta 


7.  The  environmental  settings  of  a  river,  bay,  and  delta  system 
provide  some  idealization  to  analyze  the  influencing  factors.  In  our 
study,  four  main  features  are  apparent: 

a.  The  large  input  of  fresh  water  and  sediments  from  the 
Lower  Atchafalaya  River  is  undoubtedly  the  dominant 
forcing  function  in  shaping  the  Atchafalaya  River 
Delta . 

b.  The  Atchafalaya  Delta  is  building  into  a  shallow  bay,  in 
contrast  to  the  continential  shelf  location  of  the 
Mississippi  Modern  Balize  Delta  (modern  birdfoot  delta) . 

c.  The  Atchafalaya  Bay  domain  is  constrained  within  an  area 
of  233  square  miles*  (33  miles  wide,  and  8  miles  long). 
The  average  depth  in  the  bay  is  about  5  ft,  and  the 
water  volume  is  about  3.25  x  10^®  ft^  (McAnally  and 
Heltzel  1978). 

d.  The  average  salinity  of  the  waters  in  Atchafalaya  Bay  is 
0.37  ppt  (US  Fisheries  and  Wildlife  Service  1976). 

8.  From  above  information,  we  infer  that  the  main  forces  that 
affect  the  behavior  of  Atchafalaya  River  discharge  and  subsequent 


*A  table  of  factors  for  converting  non-SI  units  of  measurement  to 
SI  (metric)  units  is  presented  on  page  7. 


delta  formation  in  the  bay  are  primarily  inertial  force  and  bottom 
friction.  Therefore  our  analytical  approach,  based  upon  the  turbu 
lent  plane  jet  theory,  to  quantify  the  delta  development  in  early 
stages  is  logical.  Our  scheme  for  characterization  of  fresh  river 
water  discharging  into  a  quiescent  bay  is  defined  in  Figure  3. 
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Figure  3.  Schematization  of  fresh  river  water  discharging 

into  a  quiescent  bay 


PART  II:  LITERATURE  REVIEW  OF  TURBULENT  JETS 


Review  of  General  Basics  of  Turbulent  Jets 

9.  Fresh  water  discharging  from  a  river  into  coastal  waters 
forms  a  turbulent  jet.  River  discharge  at  the  mouth  transports 
sediments  into  the  bay.  Mass  transport  of  these  materials  determines 
the  ultimate  distribution  of  sediments  and  the  bathymetric  changes 
near  river  mouths.  Essential  features  of  river  effluents  have  been 
summarized  briefly  in  Part  I  and  have  been  documented  thoroughly 

by  Coleman  (1976). 

10.  There  are  many  articles  and  technical  papers  on  turbulent 
jets.  The  literature  reviewed  in  this  chapter  is  limited  primarily  to 
surface  jets  and  plumes  with  bottom  friction  and  lateral  entrainment. 
Types  of  turbulent  jets 

11.  A  general  view  of  the  basics  of  turbulent  jets  is  given  by 
Pai  (1954),  Townsend  (1956),  and  Schlichting  (1968).  Much  more 
detailed  analysis  is  presented  in  the  classical  work  of  Abramovich 
(1964)  and  in  the  relatively  recent  book  by  Rajaratnam  (1976). 

12.  Turbulent  jets  are  a  special  category  of  turbulent  shear 
flows.  Depending  on  their  dominant  driving  force,  they  are  distin¬ 
guished  as  momentum  jets  or  as  plumes.  For  the  former  momentum  is  the 
predominant  factor,  while  for  the  latter  buoyancy  is  the  governing 
force.  Jets  are  also  classified  according  to  their  geometrical  shape 
as  plane  or  axisymmetrical . 

13.  There  is  usually  not  a  distinct  separation  line  between  the 
two  categories;  a  jet  is  then  of  the  mixed-type.  In  this  case,  at  the 
area  close  to  the  outlet  the  jet  is  influenced  by  the  initial  momentum 
of  the  fluid  and  is  treated  as  a  momentum  jet,  while  at  some  distance 
from  the  outlet  the  buoyant  effects  become  more  important  so  that  the 
jet  is  treated  like  a  plume. 

Surface  heated  jets 

14.  Considerable  research  in  dealing  with  surface  heated  jets 
has  been  done.  Hayashi  and  Shuto  (1967)  were  among  the  first  who 
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presented  an  approximate  theory  for  the  solution  of  thermal  jets 
discharged  horizontally  at  the  water  surface.  Many  assumptions  were 
incorporated  in  the  analysis.  Their  solution  is  valid  where  the 
Richardson's  number  is  close  to  unity. 

15.  Another  model  for  surface  thermal  jets  was  developed  by 
Hoopes,  Zeller,  and  Rohlich  (1967),  in  which  the  wind  shear  stresses 
on  the  jet  surface  and  the  entrainment  due  to  the  wind  were  in¬ 
cluded.  Their  model  is  two-dimensional  with  constant  jet  depth  and 
no  vertical  entrainment.  They  assumed  that  the  jet  spreads  linearly; 
buoyant  effects  and  current  drag  forces  are  neglected. 

16.  Motz  and  Benedict  (1970)  studied  the  problem  of  heated 
surface  jet  discharges  into  rivers.  Their  model  is  also  two-dimen¬ 
sional  with  constant  jet  thickness,  but  considers  both  vertical 
entrainment  and  drag  forces.  Although  buoyancy  was  assumed  to  induce 
vertical  entrainment,  it  was  neglected  with  respect  to  jet  spreading. 

17.  In  the  study  of  the  discharge  of  heated  water  into  deep 
receiving  waters,  Koh  and  Fan  (1970)  were  first  to  introduce  inter¬ 
facial  shear  stress  in  the  formulation  of  turbulent  jet  problems. 
Internal  hydraulic  jump  due  to  the  transition  of  the  supercritical 
issuing  jet  in  the  near-field  zone  into  a  subcritical  flow  in  the 
far-field  zone  was  discussed  in  their  study. 

18.  Shirazi  and  Davis  (1974)  have  developed  a  model  for  buoyant 
surface  jets,  which  in  essence  is  similar  to  the  one  given  by 
Stolzenbach  and  Harleman  (1971),  but  differs  in  the  fact  that  they 
used  Gaussian  similarity  profiles.  Shirazi  and  Davis  estimated  the 
coefficients  of  entrainment,  turbulent  exchange,  drag, and  shear 
through  calibration  of  field  and  experimental  data.  This  approach 
seems  less  desirable  since  many  errors  might  be  lumped  into  the 
coefficients  (Jirka,  Abraham,  and  Harleman  1975). 

Ebb  tidal  jets 

19.  Tidal  inlets  act  as  an  interface  between  estuarine  and 
coastal  waters.  Tidal  currents  near  inlets  and  estuary  mouths  play 
important  roles  in  transporting  pollutants  and  sediments.  Patterns  of 
tidal  flow  change  with  time  during  a  tidal  period.  During  ebb,  the 
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flow  on  the  ocean  side  separates  from  boundaries,  as  opposed  to  the 
nonseparated  flow  during  flood  (Ozsoy  1977).  Therefore  a  turbulent 
jet  is  often  formed  during  ebb  flows. 

20.  An  attempt  was  made  by  French  (I960)  to  model  the  ocean 
flow  patterns  during  a  tidal  period.  In  his  study,  constant  depth  and 
negligible  bottom  friction  were  assumed.  The  results  thus  obtained 
did  not  simulate  the  actual  conditions.  In  reality,  the  bottom  slope 
and  bottom  friction  may  become  important,  especially  during  ebbing 
tide  when  a  jet  flow  is  found. 

21.  Ozsoy  (1977)  incorporated  variable  bottom  topography,  bed 
resistance,  and  lateral  entrainment  into  the  jet  flow  phenomenon. 

He  described  simulated  turbulent  jet  characteristics  extensively. 

His  results  of  the  flow  patterns  at  the  vicinity  of  tidal  inlets  have 
been  compared  with  a  small  physical  model  and  good  agreement  was 
found . 

22.  Sill,  Fisher,  and  Whiteside  (1981)  investigated  deposition 
in  an  inlet  where  the  hydrodynamics  were  simulated  by  a  one¬ 
dimensional  jet,  considering  only  ebbing  tide.  Their  conclusions 
were  that  the  dimensions  of  the  equilibrium  horseshoe-shaped  shoal 
were  proportional  to  the  inlet  velocity  and  that  the  two-dimensional 
theoretical  isopachs  do  not  adequately  predict  the  shape  of  the  shoal. 
Freshwater  effluent  plumes 

23.  Bates  (1953)  suggested  that  at  most  natural  river  mouths 
freshwater  effluent  diffusion  can  be  based  upon  the  theory  of  turbu¬ 
lent  jets.  A  jet  boundary  occurs  near  the  river  mouth  as  fresh  water 
discharges  into  a  quiescent  bay.  Due  to  the  discontinuity  in  the 
velocity  of  flow,  a  zone  of  turbulent  mixing  is  established. 

2U.  Wright  and  Coleman  (1971)  suggested  that  freshwater  flow 
from  a  river  mouth,  its  deceleration,  and  consequent  sediment  deposi¬ 
tion  reflect  varying  conditions  of  outflow  inertia  and  associated 
turbulence,  bottom  fricition,  buoyancy  induced  by  density  differences, 
and  the  winds,  tides,  and  currents  of  the  receiving  basin. 

25.  Based  on  extensive  observations  and  representative  field 
data,  Wright  and  Coleman  (1971)  found  that  the  jet  expansion  rate  in 
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deep  receiving  water  can  be  expressed  by  an  explicit  function  of  the 
density  ratio  between  river  water  and  seawater  and  the  densimetric 
Froude  number.  They  concluded  that  the  flow  deceleration  is  due 
mainly  to  vertical  rather  than  lateral  mixing. 

26.  A  further  investigation  of  river  effluent  dynamics  for  the 
Mississippi  River  was  conducted  by  Wright  and  Coleman  in  1974.  In 
their  study,  they  divided  the  mouth  of  the  river  into  four  semi¬ 
discrete  regions  with  specific  morphologic  and  sedimentary  character¬ 
istics.  The  relative  contributions  of  outflow  inertia,  buoyancy, 
bottom  friction, and  marine  hydrodynamics  to  the  evolution  of  the  delta 
of  a  stratified  river  were  documented. 

27.  A  dual  treatment  of  the  effluent  plume,  both  theoretical 
and  experimental,  was  presented  by  Bowman  and  Iverson  (1978).  They 
concluded  that  the  plume  is  independent  of  the  bathymetry  of  the 
region  and  its  driving  mechanism  is  the  horizontal  pressure  gradient 
due  to  the  sloping  interface  between  the  plume  and  the  ambient  water. 
Also,  they  differentiate  plumes  from  saline  wedges  on  the  basis  that 
the  former  tends  to  entrain  and  mix  downward,  while  the  saline  wedge 
tends  to  entrain  and  mix  upward. 

Review  of  Previous  Efforts  in  Analyzing  Turbulent  Jets 


28.  The  governing  equations  of  fluid  hydrodynamics  and  substance 
conservation  constitute  the  basis  for  the  mathematical  model  for 
turbulent  jets.  The  formulation  is  a  system  of  partial  differential 
equations  of  three  to  five  equations,  that  is,  one  for  fluid  mass 
continuity;  one  to  three,  depending  on  the  nature  of  the  problem,  for 
the  momentum  components;  and  one  for  the  conservation  of  the  substance 
under  consideration. 

29.  Because  of  the  nonlinear  nature  of  the  mathematical  problem 
formulated,  no  analytical  solution  for  the  complete  system  has  yet 
been  found.  The  approximate  solution  to  this  mathematical  model  can 
be  found  by  means  of  a  numerical  technique  (finite  difference  and 
finite  element  numerical  methods)  using  a  digital  computer.  However, 


by  making  certain  simplifying  assumptions  the  system  of  partial 
differential  equations  can  be  reduced  into  a  system  of  ordinary 
differential  equations  which  can  then  be  solved  analytically  or 
numerically.  This  kind  of  approach  is  called  "integral  methods."  In 
essence,  the  integral  equations  of  mass,  momentum, and  energy  are 
utilized . 

Assumption  of  self-similarity 

30.  To  obtain  a  closed-form  solution  of  a  set  of  ordinary 
differential  equations,  the  self-similarity  of  velocity  profiles  along 
the  longitudinal  distance  is  assumed.  The  similarity  hypothesis  has  a 
firm  basis  for  the  classical  turbulent  jets  as  demonstrated  by  theory 
(Abramovich  1964;  Schlichting  1968). 

31.  Various  functional  forms  of  velocity  profiles  and  sediment 
concentration,  such  as  Gaussian  probability  distribution  (Shirazi  and 
Davis  1974),  have  been  developed  based  on  different  sets  of  hypotheses. 
Experimentally,  the  similarity  is  well  established  for  the  case  of 

free  jets.  However,  in  cases  of  attached  jets  with  interface  friction, 
the  similarity  is  under  question.  Another  approach  is  the  split  of 
the  spatial  solution  field  into  a  near-field  model  (close  to  the 
outlet)  and  a  far-field  model  (far  from  the  outlet). 

32.  Recent  experiments  by  Safaie  (1979)  suggest  that  the  simi¬ 
larity  function  depends  not  only  on  the  width  of  the  outlet  but  also 
on  the  bed  slope  and  the  aspect  ratio  (the  ratio  of  the  width  and 
depth  of  an  outlet).  Such  a  dependence  should  be  expected,  but  could 
increase  the  complexity  of  the  problem  for  reaching  an  analytical 
solution. 

Assumption  of  entrainment  velocity 

33.  A  positive  step  toward  the  mathematical  modeling  of  turbu¬ 
lent  jets  is  found  in  the  pioneering  work  of  Ellison  and  Turner 
(1959)  in  which  they  assumed  that  the  entrainment  is  proportional  to 
the  velocity  multiplied  by  an  empirical  entrainment  parameter  which  is 
a  function  of  the  Richardson  number.  Their  work  was  justified  by 
their  laboratory  experiments  for  surface  jets  and  inclined  plumes. 


34.  Hirst  (1971)  introduced  an  entrainment  function  which  was 
related  to  the  buoyancy  and  the  jet  orientation.  His  analysis  was 
based  on  the  integral  forms  of  the  conservation  of  mass,  energy,  and 
two-component  momentum  equations.  Gaussian  similarity  profiles  were 
used.  The  results  of  his  predictions  compared  with  experimental  data 
were  proven  to  coincide  satisfactorily.  The  receiving  water  was  taken 
as  quiescent  and  stratified. 

35.  A  further  analysis  of  the  entrainment  mechanism  was  con¬ 
ducted  by  Price  (1979).  Using  the  experimental  data  of  Kato  and 
Phillips  (1969)  and  of  Kantha,  Phillips,  and  Azad  (1979)  for  stratified 
and  nonstratif ied  receiving  water.  Price  computed  the  entrainment 
function  from  the  mean  buoyancy  and  momentum  equations.  In  the 
momentum  equation,  a  sidewall  friction  term  was  included  in  order  to 
incorporate  the  effects  of  the  experimental  tank  walls.  This  investi¬ 
gation  covered  a  wide  range  of  Richardson's  numbers. 

Presence  of  crosscurrents 

36.  Keffer  and  Baines  (1963)  experimentally  investigated  the 
case  of  an  axisymmetrical  turbulent  jet  subjected  to  a  crosswind. 

Their  results  concluded  that  the  position  of  the  jet  in  space  can  be 
described  by  a  single  function  of  the  entrainment  parameter  and  the 
momentum  of  both  the  jet  and  the  wind.  They  also  showed  that  the 
similarity  assumption  is  still  valid  for  the  case  of  a  crosswind. 

37.  A  comprehensive  turbulent  jet  integral  model  was  developed 
by  Stolzenbach  and  Harleman  (1971).  The  model  considered  the  cases  of 
both  nonbuoyant  and  buoyant  jets  as  well  as  buoyant  jets  in  cross- 
flows.  The  main  characteristic  of  their  analysis  was  the  separation 
of  the  jet  hydrodynamic  field  into  four  separate  regions  depending  on 
the  shear  pattern  of  the  flow.  Thus  for  the  near-field  zone,  the 
governing  equations  were  written  for  each  individual  region  separately 
and  were  then  linked  together  through  transfer  equations.  For  the  far 
field  zone, there  was  only  one  region.  By  scale  analysis  the  original 
system  was  reduced  into  a  simpler  one,  which  was  then  transformed  into 
an  ordinary  differential  system  by  utilizing  the  similarity  profiles 


in  a  polynomial  form.  The  final  system  was  solved  numerically  by  a 
fourth-order  Runge-Kutta  integration  algorithm. 

38.  Jirka,  Adams,  and  Stolzenbach  (1981)  gave  a  general  presenta¬ 
tion  of  the  buoyant  surface  jets  theory.  Based  on  dimensional  analysis 
and  some  physical  arguments,  they  defined  the  flow  using  a  set  of  in¬ 
dependent  variables,  including  the  kinematic  buoyant  flux,  the  volume 
flux,  and  a  characteristic  source  length.  Their  analysis  covers  the 
near  field  of  buoyant  jets  for  deep  or  shallow  receiving  waters.  The 
case  of  crossflow  was  included  also. 

39.  The  diffusion  of  axisymmetric  jets  into  inflowing  streams 
was  recently  investigated  by  Rajaratnam  and  Stalker  (1982).  In  their 
experiment,  the  velocity  of  the  jet  ranged  from  2  to  30  times  the 
stream  current;  experimental  results  showed  the  similarity  of  the 
velocity  profiles  except  within  the  boundary  layer  portion.  The 
various  jet  characteristics  were  correlated  to  the  excess  momentum 
thickness . 

Closed -form  analytical  solutions 

40.  Hayashi  and  Shuto  (1967)  were  the  first  to  present  an 
analytical  solution  for  the  surface  heated  discharge  problem.  In 
their  formulation,  the  following  assumptions  were  used: 

a.  In  the  momemtum  equation  the  horizontal  diffusion  terms 
balanced  the  pressure  gradient. 

b.  Similarity  existed  for  the  velocity  profiles  in  hori¬ 
zontal  and  vertical  planes. 

c.  The  entrainment  rate  was  proportional  to  a  charac¬ 
teristic  velocity. 

d.  The  turbulent  diffusion  coefficients  were  constant. 
Furthermore,  they  neglected  the  advective  terms,  the  vertical  veloc¬ 
ities,  and  the  shearing  forces  at  the  surface  and  the  bottom.  A  closed 
form  solution  was  obtained  for  the  velocities  by  a  biharmonic  stream 
function  equation  where  the  entrainment  was  taken  as  zero.  Their 
solution  holds  true  where  the  Richardson  number  is  close  to  unity. 

41.  An  analysis  of  the  surface  heated  jets  at  small  Richardson's 
numbers  has  been  done  by  Engelund  and  Pederson  (1973).  Their  main 
assumptions  were: 
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a.  Similarity  profiles  for  velocities  and  densities. 

b.  Linear  density  variation  along  the  depth. 

c.  Hydrostatic  pressure. 

d.  Shear  stresses  existing  only  in  the  horizontal  planes. 

By  considering  that  the  longitudinal  momentum  dominates  the  pressure 
forces,  they  derived  a  system  of  ordinary  equations  from  which  a 
closed-form  solution  was  obtained. 

42.  In  1976  Engelund  improved  this  earlier  model.  Based  on 
their  original  system  of  equations,  he  gave  a  second  order  closed-form 
solution  for  the  near  field  and  moderate  Richardson's  numbers,  using  a 
perturbation  technique.  His  solution  is  the  only  analytical  one  which 
does  not  give  similar  profiles. 

43.  Abraham  (1976)  presented  an  analytical  form  of  the  axisym- 
metric  momentum  jets  and  plumes  in  stagnant  and  flowing  receiving 
waters.  He  described  the  limits  of  the  jet  diffusion  theory  based  on 
the  similarity  assumption  and  the  entrainment  concept.  A  compre¬ 
hensive  study  was  done  by  Policastro  and  Dunn  (1976)  on  the  integral 
models  of  surface  thermal  plumes.  Their  investigation  is  thorough  and 
outlines  the  advantages  and  disadvantages  of  the  various  models. 

44.  A  model  for  nonbuoyant  jets  in  shallow  receiving  waters  for 
the  case  of  sediment  transport  was  developed  by  Ozsoy  (1977).  He 
integrated  the  shallow-water  wave  equations  along  the  jet  width, 
including  the  lateral  entrainment  and  bottom  friction.  Using  the 
similarity  functions  for  near-  and  far-field  zones  as  given  by 
Stolzenbach  and  Harleman  (1971),  Ozsoy  obtained  analytical  solutions 
to  the  jet  equations. 

45.  The  articles  discussed  briefly  in  this  chapter  suffice  only 
for  a  general  review  of  turbulent  jets.  It  is  not  an  exhaustive  list. 
Our  attempt  is  to  cover  the  evolving  theory  of  turbulent  jets  and  the 
developing  stage  of  analysis  in  jet  phenomenon. 

46.  Our  main  task  is  to  formulate  the  freshwater  discharge  into 
a  quiescent  bay  as  a  two-dimensional  plane  jet,  to  derive  an 
analytical  solution  for  the  governing  equations  of  fluid  dynamics  and 


PART  III:  DEVELOPMENT  OF  THE  ATCHAFALAYA  RIVER  DELTA 


General  Description 

47.  The  development  of  the  river  delta  in  Atchafalaya  Bay  since 
1950  has  provided  the  opportunity  to  document  the  evolution  of  two  new 
Mississippi  delta  lobes,  the  Lower  Atchafalaya  River  Delta  and  the  Wax 
Lake  Delta.  Numerous  descriptive  studies  (Garrett,  Hawxhurst,  and  Miller 
1969;  Shlemon  1972;  Cratsley  1975;  Roberts,  Adams,  and  Cunningham  1980; 
Adams  and  Baumann  1980;  Van  Heerden  1980;  and  Van  Heerden,  Wells,  and 
Roberts  1981)  have  been  conducted  in  the  bay,  which  form  the  basic 
foundation  for  this  research. 

48.  Delta  development  is  primarily  the  product  of  an  interplay 
between  river  sediment  input  and  reworking  by  physical  processes  in 
the  receiving  basin  (Wright  and  Coleman  1974).  However,  the 
Atchafalaya  Delta  is  fundamentally  different  from  the  Mississippi 
Balize  Delta  (modern  bird-foot  delta).  The  Atchafalaya  Delta  is 
building  into  a  nonstratif ied  shallow  bay  protected  by  a  series  of 
discontinuous  oyster  shell  reefs,  as  shown  in  Figure  4  (Shlemon 
1972).  This  reef  chain,  known  as  the  Point  Au  Fer  Reef,  forms  the 
Atchafalaya  Delta's  seaward  margin,  in  contrast  to  the  continental 
shelf  location  of  the  Mississippi  modern  bird-foot  delta  (Van  Heerden 
1980).  The  Point  Au  Fer  Reef,  including  its  submarine  extension,  was 
about  10  miles  long  before  it  died  out  in  the  late  1960's,  due  to  the 
increasing  influx  of  fresh  water  and  sediment  into  Atchafalaya  Bay 
(Shlemon  1972). 

49.  Cratsley  (1975)  showed  that  the  quantity  and  size  distribu¬ 
tion  of  sediment  available  to  the  Atchafalaya  Delta  are  directly 
related  to  the  modern  history  of  the  Atchafalaya  Basin  and  River.  The 
Atchafalaya  system  is  presently  river-dominated;  average  salinity  of 
the  waters  in  Atchafalaya  Bay  is  0.37  ppt  (US  Fisheries  and  Wildlife 
Services  1976).  Salt-wedge  intrusion  does  not  appear  to  signifi¬ 
cantly  affect  the  system. 


Figure  4.  Location  map  of  the  Lower  Atchafalaya  area  (Shlemon  1972) 


50.  This  study,  as  outlined  in  PART  I,  is  focused  on  use  of  a 
freshwater  discharge  into  a  quiescent  bay  for  the  prediction  of  the 
extent  of  area  and  mass  of  delta  growth  in  early  stages.  The  evolu¬ 
tionary  history  of  the  Atchafalaya  Delta,  gathered  from  literature, 
serves  as  the  basis  of  our  evaluation  of  the  phenomena  to  be  analyzed. 
The  pertinent  information  on  river  discharge,  sediment  characteristics, 
and  bay  bathymetry  are  reviewed  and  arranged  in  chronological  order  in 
the  following  sections. 


Historical  Development 


Prior  to  1950' s 

51.  The  Atchafalaya  River  system  flowed  through  a  broad  basin 
characterized  by  extensive  freshwater  swamps  and  numerous  small  lakes. 
Prior  to  the  early  1950' s,  most  of  the  sediments  were  trapped  in  the 
catchment  basin  before  they  reached  Atchafalaya  Bay.  The  bottom  con¬ 
figuration  of  the  bay  virtually  remained  unchanged.  The  bay  depth  was 
maintained  at  a  constant  depth  of  6  ft  (Shlemon  1972).  Very  little 
sediment  was  deposited  in  the  bay.  Prodelta  clays  and  silty  clays 
were  accreted  on  the  continental  shelf  off  the  Atchafalaya  Bay 
(Cratsley  1975). 

1952  to  1962 

52.  From  1952  to  1962,  as  the  diversion  of  Mississippi  River 
flow  through  the  Atchafalaya  River  system  increased  steadily  and  as 
the  catchment  basins  were  filled,  accelerated  sedimentation  in 
Atchafalaya  Bay  marked  the  beginning  of  a  subaqueous  delta  (Cratsley 
1975).  As  displayed  in  Figure  5,  Grand  Lake  and  Six  Mile  Lake,  which 
had  been  natural  settling  basins  for  coarse-grained  sediments  (silt 
and  sand),  were  rapidly  filled  by  deltaic  deposits.  These  catchment 
basins  are  now  directly  routing  sediment  through  the  natural  Lower 
Atchafalaya  River  Outlet  and  the  man-made  Wax  Lake  Outlet. 

53.  Shlemon  (1972)  indicated  that  about  47  square  miles  of 
the  bay  bottom  had  been  covered  by  at  least  a  1-ft  thickness  of 
sediment  by  1962;  over  6  ft  of  fill  occurred  just  south  of  Shell 
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Figure  5.  Lacustrine  delta  accretion  in  the  Lower  Atchafalaya 
Basin,  1917-1960  (Shlemon  1972) 
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Island,  as  shown  in  Figure  6  and  Table  1.  In  Figure  6,  the  delta 
front  environment  in  the  Wax  Lake  Outlet  mouth  area  is  much  less 
extensive  than  that  of  the  Lower  Atchafalaya  River  mouth  area,  re¬ 
flecting  the  relative  discharges  from  the  two  outlets  during  the 
period  of  1952  to  1962  (Cratsley  1975). 

1965  to  1967 

54.  In  the  late  1960's,  the  increasing  discharge  of  the 
Atchafalaya  River  and  the  increasing  rate  of  suspended  sediment 
transport  were  responsible  for  the  formation  of  the  distributary  bars 
at  the  mouths  of  the  Wax  Lake  Outlet  and  the  Lower  Atchafalaya  River 
(Cratsley  1975).  The  average  monthly  discharge  at  the  latitude  of 
the  outlets  during  the  period  1965-1967  and  the  corresponding  sus¬ 
pended  load  transported  through  the  outlets  were  summarized  by  Garrett, 
Hawxhurst,  and  Miller  (1969),  and  Cratsley  (1975)  as  reproduced  in 
Figure  7.  The  combined  average  annual  freshwater  flow  through  the 
outlets  was  about  165,000  cfs.  High  riverflows  occurred  from 
January  through  June,  reaching  a  maximum  of  325,000  cfs  in  May  and  a 
minimum  of  73,000  cfs  in  September.  The  textural  composition  of  the 
suspended  load  was  25  percent  sand  and  75  percent  silt  and  clay 
(Garrett,  Hawxhurst,  and  Miller  1969). 

1967  to  1972 

55.  The  period  1967  to  1972  was  characterized  by  the  expansion 
of  the  delta  front  environment  throughout  Atchafalaya  Bay,  the  estab¬ 
lishment  of  distributary  mouth  bars  in  the  bay,  and  the  rapid  prograda¬ 
tion  of  the  distal  bar  (Cratsley  1975).  Very  definite  subaerial 
delta  lobes  appeared  in  1972.  These  initial  subaerial  exposures  were 
shoals,  composed  largely  of  sediment,  extending  from  the  Atchafalaya 
River  Outlet  to  Point  Au  Fer  Shell  Reef. 

56.  The  rapid  expansion  of  the  delta  front  environment  in  the 
Wax  Lake  Outlet  channel  mouth  area  is  shown  in  Figure  8,  the  1972 
bathymetric  map  prepared  by  the  US  Army  Engineer  District,  New  Orleans 
(Adams  and  Baumann  1980).  A  more  detailed  evaluation  of  the  prograda¬ 
tion  of  the  delta  front  environment  was  made  by  Cratsley  (1975)  by 


Figure  6.  Isopach  map  showing  area  and  thickness  of  deltaic  fill  in 
Atchafalaya  Bay,  1952-1962  (Shlemon  1972) 


measuring  the  changes  in  bottom  topography  along  the  six  profile  lines 
as  delineated  in  Figure  9  and  plotted  in  Figure  10. 

57.  The  annual  flow  in  the  Atchafalaya  River  at  Siininesport , 
Louisiana  (near  the  diversion  point  in  the  upper  basin),  from  1967  to 
1972,  was  208,000  cfs,  and  the  average  annual  peak  flow  was  303,000 
cfs  (Adams  and  Baumann  1980).  Adams  and  Baumann  (1980)  indicated 
that  the  increase  in  discharge  during  the  period  of  1967  to  1972 
occurred  predominantly  during  traditionally  low  flow  months;  the  peak 
discharges  during  the  spring  months  were  not  great  enough  to  increase 
sediment  load  entering  the  bay.  The  average  annual  sediment  load 
delivered  to  Atchafalaya  Bay  was  about  63  x  10®  tons  for  the  period 
1965-1971  (USACOE  1974). 

1973  to  1975 

58.  The  1973-1975  years  were  abnormally  high-water  years 
compared  with  the  past  20  years  of  flow  records  at  Simmesport, 
Louisiana,  as  displayed  in  Figure  11  (Van  Heerden  1980).  During 
1973-1975,  flows  averaged  315,000  cfs  at  Simmesport;  peak  flows  of 
over  700,000  cfs  occurred  in  April  1973  and  over  600,000  cfs  in  April 
1975. 

59.  Similar  high-flow  averages  and  peak  flow  were  recorded  at 
Morgan  City,  on  the  Lower  Atchafalaya  River,  during  1973-1975.  Both 
the  discharge  and  suspended  load  at  Morgan  City  are  displayed  in 
Figure  12  (Roberts,  Adams,  and  Cunningham  1980).  Peak  flows  at  Morgan 
City  of  over  600,000  cfs  occurred  in  May  1973,  and  the  normal  300,000- 
cfs  peak  flows  were  exceeded  during  8  months  of  1973-1975. 

60.  Accompanying  these  abnormally  high  discharges  were  un¬ 
usually  high  concentrations  of  sediment  carried  as  suspended  load 
(25  percent  sand,  75  percent  silt  and  clay).  The  annual  suspended 
sediment  load  reaching  Atchafalaya  Bay  during  the  three  high-water 
years  was  about  123  x  10®  tons.  The  sediment  budget  and  size 
characteristics  during  the  periods  of  1965-1971  and  1973-1975  are 
listed  in  Table  2  (Roberts,  Adams,  and  Cunningham  1980). 


Figure  9.  Index  map  for  bottom  topography  profile  lines  (Cratsley  1975) 


Discharg* 


Figure  12.  Discharge  and  suspended  load  at  Morgan  City  near  Lower  Atchafalaya  River  Outlet 

1973-1975  (Roberts,  Adams,  and  Cunningham  1980) 


1976  to  1978 


61.  Although  the  years  1976  through  1978  were  considered  as 
normal  flood  years  (Figure  11),  the  previous  three  abnormally  high- 
water  years  had  played  an  important  role  in  the  rapid  development  of 
the  subaerial  delta  phase.  During  1976-1978,  the  distributary  mouth 
bar  extended  seaward  and  evolved  into  a  complex  branching  network 
characteristic  of  deltas  where  river  mouths  are  frictionally  dominated 
and  are  gradually  building  into  low-energy,  shallow-water  environments 
(Wright  and  Coleman  1974). 

62.  Bathymetric  data  taken  in  1977  by  USACOE  and  adjusted  to 
the  1975  msl  datum  (Adams  and  Baumman  1980)  indicated  that  an  esti¬ 
mated  6.55  square  miles  of  new  land  has  developed  above  msl  (Figure  13). 
Above  the  -1  ft  datum,  which  represents  the  mean  low  tide  level,  a 
calculated  15.8  square  miles  of  new  subaerial  land  with  an  approximate 
width  of  6.8  miles  had  been  added  to  Atchafalaya  Bay  over  the  period 
1967-1977  (Roberts,  Adams,  and  Cunningham  1980). 

1979 

63.  Another  major  flood  occurred  in  1979.  A  peak  flow  of  over 
500,000  cfs  was  recorded  in  April  1979  (Figure  11).  Roberts,  Adams, 
and  Cunningham  (1980)  concluded  that  suspended  sediment  transport 
during  floods  was  responsible  for  the  abrupt  increases  in  subaerial 
delta  growth. 

64.  Through  using  satellite  imagery,  color  infrared  photog¬ 
raphy,  and  digital  current  meter  data,  Wells  and  Kemp  (1981)  provided 
estimates  on  the  suspended  sediment  concentrations  within  Atchafalaya 
Bay.  These  average  about  300  mg/1  and  range  from  250  to  400  mg/1. 

Long-Term  Future  Projection 

1970  to  2020 

65.  Shlemon  (1972),  using  sediment  load  measurements  obtained 
in  the  outlets,  outlined  the  probable  future  configuration  of  the 
Atchafalaya  Delta  by  the  year  2020  (Figure  14).  An  average  growth 
rate  was  inferred  and  is  plotted  in  Figure  15.  A  straight  line 
projection  of  bay  filling  until  the  year  2020  will  yield  an  estimated 
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Inferred  50-year  (1970-2020)  average  growth  rate  of  the  Atchafalaya 
River  Delta  if  unmodified  by  man 


filling  rate  of  about  7  square  miles  per  year,  covering  an  area  of  350 
square  miles.  Assuming  that  the  peak  of  subaerial  growth  will  be 
reached  by  1990,  a  conservative  growth  rate  will  be  5.5  square  miles 
per  year,  with  resultant  production  of  nearly  290  square  miles  of  new 
land  by  the  year  2020. 

66.  Roberts,  Adams,  and  Cunningham  (1980)  pointed  out  that  unless 
an  abnormal  number  of  peak  floods  such  as  those  of  the  1973-1975  period 
occur  during  the  next  two  decades,  Atchafalaya  Bay  probably  will  not  be 
filled  until  after  the  turn  of  the  century.  They  further  estimated 
that  the  sand-dominated  phase  of  the  delta  will  probably  cover  an  area 
of  over  50  square  miles  before  shifting  its  locus  of  deposition  to  the 
shelf  seaward  of  Che  Point  Au  Per  Shell  Reef. 

67.  Van  Heerden,  Wells,  and  Roberts  (1981)  projected  that  the 
Atchafalaya  Delta  should  prograde  more  rapidly,  form  thin  sand  bodies, 
and  eventually  cover  a  wider  area,  much  like  the  Lafourche,  St.  Bernard, 
and  Teche  delta  lobes. 

1977  to  2027 

68.  A  statistical  approach  to  predict  the  future  growth  of  the 
Atchafalaya  River  delta,  based  on  historical  deposition  trends  in 
Atchafalaya  Bay,  has  been  presented  by  Letter  (1982).  He  developed  a 
regression  model  that  correlates  deposition  rates  with  river  dis¬ 
charge,  sediment  yield,  water  depth,  and  the  delta  mass  centroid. 

Using  the  1977  bathymetry  as  an  initial  condition,  the  model  is 
applied  to  a  50-year  hydrograph  at  10-year  increments  of  prediction. 

The  results  of  the  regression  model  showed  that  within  50  years  the 
delta  will  evolve  gulfward  of  Eugene  Island,  the  gulfward  limit  of  the 
bay. 

69.  Figure  16  shows  the  predicted  condition  of  the  Atchafalaya 
Bay  in  the  year  2027.  The  total  volume  of  the  deposited  sediments  is 
estimated  at  58  billion  ft^,  and  the  delta  mass  volume,  based  on  -3  ft 
NGVD  (National  Geodetic  Vertical  Datum),  is  about  17.6  billion  ft^ 
(Letter  1982). 
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Estimated  delta  mass  and  deposited  sediment  volume 
for  the  Atchafalaya  Delta,  1970-2020 


CUMULATIVE  VOLUME  DEPOSITED  IN  BILLIONS  OF  CUBIC  FEET 


1980  to  2030 


70.  Long-term  predictions  based  on  the  morphologic  development 
of  the  Atchafalaya  River  delta  and  generic  analysis  of  existing  deltas 
have  been  completed  recently  by  Wells,  Chinburg,  and  Coleman  (1984). 
Utilizing  maps,  charts,  aerial  photographs,  and  LANDSAT  imagery,  they 
have  examined  10  similar  deltas  and  subdeltas  worldwide,  and  have 
projected  the  rate  of  growth  and  configuration  of  subaerial  land  in 
Atchafalaya  Bay  to  the  year  2030. 

71.  The  requirements  for  similar  deltas  that  closely  resemble 
the  Lower  Atchafalaya  River  delta  were  defined  by  Wells,  Chinburg,  and 
Coleman  (1984)  in  their  generic  analysis.  The  requirements  were  low 
wave  energy,  low  tidal  energy,  a  shallow  receiving  basin,  and  high 
suspended  sediment  load.  The  results  from  their  study  indicate  that 
the  subaerial  land  area  in  Atchafalaya  Bay  by  the  year  2030  will  range 
from  150  km^  (59  square  miles)  to  337  km^  (132  square  miles),  with 
208  km^  (81  square  miles)  representing  the  expected  land  area  in  50 
years  under  average  flood  conditions. 

72.  Approximately  14  x  10®  m^  (495  x  10®  ft^)  of  sediment  per 
year  is  retained  in  Atchafalaya  Bay  (Wells,  Chinburg,  and  Coleman 
1984).  Growth  prediction  curves  for  subaerial  land  in  the  bay  were 
constructed  by  Wells,  Chinburg,  and  Coleman  (1984),  as  shown  in 
Figure  17.  Figure  18  shows  the  configuration  of  land  in  the  bay  in 
the  year  2030  based  on  the  range  of  predicted  rates  of  growth. 

73.  Table  3  summarizes  the  development  of  the  Atchafalaya  River 
delta  in  chronological  order  as  gathered  from  various  literature. 
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Growth  curves  predicting  subaerial  land  in  Atchafalay 


PART  IV:  GOVERNING  EQUATIONS  OF  TURBULENT  PLANE  JETS 
IN  SHALLOW  RECEIVING  WATERS 


General  Description  of  River-Bay  Systems 

Ik.  River  mouths,  with  their  unique  features,  are  very  complex 
physical  systems.  Because  they  provide  a  natural  link  between  inland  and 
sea,  most  of  the  world's  civilizations  and  metropolitan  cities  have  been 
developed  close  to  these  areas.  Therefore  conflicting  interests  re¬ 
garding  ecologic,  economic,  recreational,  and  transportation  issues  are 
closely  interrelated  with  these  coastal  domains,  creating  a  challenge 
where  both  human  and  natural  forces  are  involved. 

75.  The  main  factors  affecting  the  processes  of  river  mouth  delta 
formation  are  river-bay  hydrodynamics,  geomorphology,  climate,  and  human 
activities  (Coleman  1976).  Due  to  the  practical  importance  of  river 
delta  development,  many  scientific  efforts  have  been  placed  on  the 
investigation  of  basic  laws  and  principles  that  govern  the  behavior  and 
response  of  fresh,  sediment-laden  river  water  as  it  issues  into  the 
receiving  salt  waters  of  the  sea. 

76.  The  overall  river-bay  system  is  a  time-dependent,  three- 
dimensionsal  phenomenon,  influenced  by  a  vast  number  of  different 
deterministic  and  stochastic  parameters.  To  develop  a  general  mathe¬ 
matical  model  would  be  a  formidable  task ;  thus  based  on  field  and 
laboratory  observations  of  similar  phenomena,  certain  assuiiiptions  are 
made  and  predominant  parameters  with  deterministic  characteristics  are 
used  to  reduce  the  problem  considerably  and  to  make  it  suitable  for 
theoretical  formulation  and  analytical  approach. 

Formulation  of  River  Discharge  into  a  Bay 

77.  Physically,  a  river-bay  system  can  be  regarded  as  a  plane 
water  jet  issuing  into  a  large  receiving  body  of  water.  Mathematically, 
the  jet  hydrodynamics  of  a  river-bay  system  can  be  expressed  by  a  system 


of  time-dependent,  nonuniform,  incompressible,  f ree-surface ,  three- 
dimensional  turbulent  flow  equations  (Schlichting  1968).  In  general, 
this  is  a  nonlinear,  hyperbolic-type  partial  differential  system  of  four 
equations,  one  for  mass  continuity  and  three  for  momentum  balance.  The 
independent  variables  are  the  three  Cartesian  axes  (x^,  x^,  x^)  and  the 
time  (t),  as  shown  in  Figure  3,  while  the  dependent  variables  are  the 
three  velocity  components  parallel  to  the  three  axes  and  the  water 
elevation  (q) .  In  many  practical  cases  this  system  can  be  reduced  into  a 
two-  or  one-dimensional  model  and  still  be  able  to  describe  the  physical 
phenomena  adequately. 

78.  The  general  form  of  the  equation  of  continuity,  or  conserva¬ 
tion  of  mass,  for  an  incompressible  flow  (Schlichting  1968)  is  written 
as 


i  =  1,  2,  3  (1) 


where  u^  =  the  velocity  components,  and  x^  =  the  Cartesian  coordinates. 

79.  The  equations  of  momentum  are  derived  from  the  Navier-Stokes 
equation  (Schlitching  1968)  and  are  expressed  as 


Du.  1  8p  1  3t., 

_ 1  _  _  _  ^  ^ 

Dt  p  8x .  p  8x, 

k 


k  =  1,  2,  3  (2) 


where  the  symbol  D/Dt  denotes  the  total  derivative,  that  is,  the  equiva¬ 
lence  of  an  operator. 
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p  =  the  density,  p  =  the  pressure,  T . .  =  the  stress  components,  and  b.  = 

1  iC  X. 

the  body  forces. 
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Shallow-Water  Hydrodynamic  Equatioas 


80.  The  general  system  can  be  reduced  into  a  simpler  form  for  the 
case  of  shallow  waters  as  pertaining  to  the  Atchafalaya  River-Bay  system. 
The  basic  assumption  is  the  hydrostatic  pressure  distribution,  and  by 
neglecting  the  vertical  stress  components  and  the  vertical  acceleration, 


the  x_  momentum  balance  in  Equation  2  reduces  to 


1  8p 


“  "  ■  P  3*3  ■  * 


(4) 


or 


p  =  pg(n  -  x^)  +  Pq 


(5) 


where  n  =  free  surface  elevation  from  the  reference  datum  and  p^  = 


atmospheric  pressure.  The  reference  datum  coincides  with  the  mean  sea 
level  (Figure  3). 

Equation  of  continuity 

81.  The  free  surface  can  be  expressed  as  =  r|(Xj,  x^,  t). 
Differentiating  this  with  respect  to  time  yields 
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Similarly,  if  x^  -  -h(Xj,  X2)  is  the  distance  from  the  reference  datum  to 
the  bottom,  then 


U3  =  -  u^ 
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Integration  of  Equation  1  with  respect  to  the  x^-axis  gives 
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-h 


-  dx- 
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C  ^“2 
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a.,  ^ 

J  ax 

-h 

Leibnitz ' 

s  rule  of  : 

6  and  7 , 

Equation  8 
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J  ax¬ 
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a^..a^  ..an„ 
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Furthermore,  due  to  the  shallowness  of  the  bay  waters,  the  velocities  are 
assumed  to  be  uniform  along  the  depth.  Then  Equation  9  can  be  written 
as 


a(h  +  n)u,  a(h  +  n)u^  an 

- -  +  - =  +  —  =  0 


(10) 


ax 


1 


ax^  at 


Equation  10  is  the  final  form  of  the  continuity  equation  in  a  two- 
dimensional,  time-dependent  shallow-water  horizontal  domain. 

Equation  of  momentum 

82.  After  introducing  the  shallow-water  approximation  and 
expanding  the  total  derivative.  Equation  2  becomes 
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(11) 
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i  =  1,  2 
k  =  1,  2,  3 

Let  the  time  average  of  the  velocity  components  be  represented  by 


L^i  =  ^2’  ^  ^2’ 


(12) 


where  =  mean  velocities  and  =  perturbation  velocities  with  zero 
mean  value.  Also  let  the  only  mass  force  be  the  Coriolis  force.  Thus, 
for  the  northern  hemisphere, 


b^  =  2  Q  sin»5  (13a) 

b^  =  -2  fi  sinfi  Uj  (13b) 


where  Q  =  angular  velocity  of  the  earth,  and  ^  -  geographical  latitude  of 
the  river  mouth  of  the  site.  The  shearing  forces  can  be  approximated  as 


(14) 


where  p  =  molecular  viscosity. 

83.  Substituting  Equation  12  for  the  perturbated  velocities. 
Equations  13a  and  13b  for  the  Coriolis  forces,  and  Equation  14  for 
the  approximated  shearing  forces  into  Equation  11  and  time -averaging 
it,  one  arrives  at 
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1.  11.  IZ.OlJ  ^  rs  ■  J  ~ 

-  +  -  +  -  +  -  -  2  \l  sinp  u.. 


1  ap  M  /  a  u  a  u  a  i 

- +  -( - 2  +  - 2“  ^ 

p  ax^  p  \  axj  ax^  Bx^‘‘ 


au^u^  au^u^ 


(15a) 


3u  3(u  u  )  3(u  u  )  (a  u  )  . 

-  +  -  +  -  + - +  2  Q  sinji  u. 


1  ap  M  /  ^^^2  ^^“2  ^^“2 

p  ax.,  p  \  ax,^  ax-^ 


^  ^  3u^u’ 


(15b) 


84.  The  last  two  terras  of  Equation  15a  and  15b  can  be  grouped 
together  as  total  stress  terras  designated  by  the  eddy  viscosity  stresses 
in  both  horizontal  and  vertical  components.  Equations  15a  and  15b 
become 


au  a(u  u  )  a(u  u  )  a(u  u  ) 

— t  +  - LJ_  +  - L_^  +  - 

at  ax,  ax-  ax„ 


2  Q  sinjJ  u.. 


p  ax . 


,  +  !Jii 

^  \  ax^^  ax^^ 


(16a) 
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-  +  -  +  -  +  -  +  2  12  sinp  u 
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axi 

1  ap 

p  ax. 


ax. 
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+  £ 


(16b) 


where  =  horizontal  eddy  viscosity  coefficient,  and  =  vertical  eddy 
viscosity  coefficient.  The  bars  (-)  were  dropped  for  convenience. 

85.  Neglecting  the  vertical  velocity  component  (u^))  substituting 
the  pressure  p  from  Equation  5,  integrating  Equation  16  from  bottom 
to  surface  along  the  x^-axis,  dividing  by  h  +  r)>  ^nd  assuming  vertically 
averaged  velocities,  Equations  16a  and  16b  become 


and 


au  a(u  u  )  a(u  u,) 

-  + - +  -  -  2  Q  sinji  u^ 

at  axj  ax2 


an 

au^ 

ax^ 

3X27 

h+n 

ax3 

(17a) 


au  a(u  u  )  a(u  u  ) 

—  +  - +  2  n  sinji  u 

at  ax^  ax^  ^ 


n 


(17b) 

where  u^  and  u^  now  stand  for  the  average  velocities  over  the  depth  in 
the  horizontal  domain. 
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86.  The  last  term  in  Equations  17a  and  17b  is  the  expression 
for  the  tangential  stresses  at  the  surface  (x^  =  rj)  and  at  the  bottom  (x^ 
=  -h) .  Experimental  studies  in  a  one-dimensional  rivertlow  resulted  in 
an  empirical  formula  for  the  bottom  stresses 


2 

u 


z 


where  =  Chezy's  coefficient  of  friction  and  u  =  the  mean  cross- 

sectional  velocity.  The  value  of  depends  on  various  geometrical  and 

flow  parameters.  It  is  suggested  that  this  coefficient  be  evaluated 

through  calibration  from  actual  field  data.  Experience  shows  that  the 

k  k 

coefficient  usually  varies  between  45  m'/sec  and  70  m*/sec. 

87.  Similarly,  the  air-sea  interaction  stresses  (t^)  can  be 
approximated  by  a  formula  like  that  of  Equation  18  and  can  be  written 
as 

T  =  p  y^w^  (19) 

s  a  * 

where  p  =  density  of  the  air,  y  =  constant  coefficient  and  w  =  wind 

3 

velocity  at  some  reference  height  from  the  water  surface.  Experimentally 

2 

it  has  been  found  that  y  is  a  function  of  the  wave  form,  with  a  value 
close  to  0.0026.  Further  remarks  on  the  bottom  and  surface  shear 
stresses  can  be  found  in  Dronkers  (1964)  and  Nihoul  (1975). 

88.  Neglecting  the  effects  of  shear  and  turbulent  mixing,  and  with 
the  introduction  of  surface  stress  (Equation  19)  and  bottom  stress 
(Equation  18),  Equations  17a  and  17b  can  be  written  as 
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where  =  the  wind  velocity  component  along  the  x^-axis,  and  w^  =  the 
wind  velocity  component  along  the  x^-axis.  Equations  20a  and  20b  are 
the  momentum  balance  equations  for  a  two-dimensional  shallow-water  wave 
hydrodynamic  model. 

89.  For  the  completeness  of  the  mathematical  model,  the  proper 
conditions  must  be  defined  at  the  spatial  domain  boundary,  and  the 
conditions  at  the  initiation  time  of  the  phenomenon  must  also  be 
provided. 


General  Description  of  River-Delta  Interaction 


90.  River  delta  development  is  primarily  the  interaction  of  river 
sediment  input  and  the  physical  processes  of  the  receiving  basin 
(Coleman  1976).  The  large  input  of  fresh  water  that  carries  sediments 
from  the  Lower  Atchafalaya  River  is  the  dominant  force  in  shaping  the 
Atchafalaya  River  delta. 

91.  The  nature  of  the  sediments  plays  an  important  role  in  the 
river  delta  system.  The  presence  of  cohesive  or  noncohesive  sediments, 
and  the  critical  shear  stress  for  erosion,  deposition,  suspension, or 


consolidation  are  controlling  features  in  delta  formation.  Cohesive 
sediments  are  more  complicated  to  deal  with,  since  they  are  controlled 
not  only  by  the  hydrodynamic  forces  but  also  by  the  electrochemical 
forces  (Krone  1978) . 


Sediment  Transport  and  Delta  Formation 


92.  Freshwater  effluent  from  river  mouths  carries  sediments  in 
suspension.  The  diffusion  of  these  materials  and  mixing  with  the  ambient 
bay  water  determine  their  ultimate  distribution.  Based  on  the  hydro- 
dynamic  aspects  presented  in  previous  sections,  the  study  of  turbulent 
jet  diffusion  processes  in  shallow  water  is  formulated  in  the  following 
sections . 

Shallow-water  mass  transport  system 

93.  When  the  distribution  of  a  physical  property  or  substance 
(sediment,  salinity,  temperature,  or  chemical  wastes)  that  is  carried  by 
the  jet  needs  to  be  studied,  an  additional  equation  is  then  utilized. 

This  is  the  convection-diffusion  equation  with  proper  source  or  sink 
terras.  The  equation  accounts  for  the  mass  conservation  and  is  a  partial 
differential  equation  of  the  parabolic  type,  having  as  a  dependent 
variable  the  concentration  of  the  substance  under  consideration. 

Equation  of  mass  conservation 

94.  In  general,  the  conservation  of  mass  of  a  substance  is  ex¬ 
pressed  by  the  convective-diffusion  equation  (Ariathurai,  MacArthur, 
and  Krone  1977)  in  the  form 

^  +  S  i  =  1,  2  (21) 

where  c  =  substance  concentration,  D.  =  the  molecular  diffusion  coeffi- 

1 

cient,  and  S  =  the  proper  source  and/or  sink  term. 

95.  Assuming  again  a  time  mean  and  a  perturbated  value  for  the 

velocities  and  concentrations,  that  is,  u.  =  u.  +  ul ,  and  c  =  c  +  c' , 
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substituting  these  terms  into  Equation  21,  and  time-averaging  it  one 
arrives  at 


3c 

3c 

3(u'.c') 

A  ^  — 
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1  u 
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-  r  u . 
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The  deviation  term  can  be  approximated  as 


(22) 


3(u:E') 

1 

3x . 

1 


(23) 


where  =  the  eddy  diffusion  coefficient. 

96.  Combining  the  expressions  for  molecular  and  eddy  diffusion, 
and  incorporating  their  coefficients  into  a  single  term  (E^),  Equation 
22  becomes 


+  S 


(24) 


where  the  bars  have  been  dropped  for  simplicity.  Equation  24  is  the 
general  equation  of  the  conservation  of  a  substance;  when  it  refers  to 
conservation  of  sediments,  the  term  S  stands  for  the  processes  of 
erosion  and/or  deposition. 


PART  V:  ANALYTICAL  TECHNIQUES  AND  NUMERICAL  SOLUTIONS 


Simplifying  Assumptions  Applying  to  the 
Atchafalaya  River-Bay  System 

97.  The  system  of  Equations  10,  20,  and  24  derived  in  PART  IV 
is  a  complicated  system  of  partial  differential  equations  that  can  be 
solved  by  means  of  a  numerical  technique  through  a  digital  computer. 
Analytical  solutions  can  be  achieved  only  in  the  case  where  the  general 
equations  are  simplified  considerably  under  certain  assumptions.  These 
assumptions  can  be  derived  from  the  specific  features  of  each  river-bay 
system  and  the  characteristic  properties  of  the  jet  itself.  However, 
misuse  of  the  simplifications  may  lead  to  erroneous  solutions  of  little 
practical  value.  Thus  an  extended  and  detailed  knowledge  of  the  physical 
system  under  consideration  is  required  so  that  the  limitations  of  the 
validity  of  the  solution  can  be  well  understood. 

98.  Regarding  the  Atchafalaya  River-Bay  system,  on  a  first 
approximation  basis,  the  following  assumptions  were  used: 

a.  Shallow  receiving  waters  and  velocities  are  uniform  over 
depth,  where  friction  is  a  predominant  factor. 

b.  Well-mixed  conditions  with  no  density  stratification. 

c.  Negligible  density  difference  between  issuing  and 
receiving  waters,  that  is,  a  nonbuoyant  jet. 

d.  Very  small  wave  height,  rj>  in  comparison  to  the  depth, 
h,  so  that  n  ~  0. 

e.  Negligible  Coriolis  forces  and  wind  stress  effects. 

f.  Bell-shaped  similarity  profiles  for  the  velocities 
and  sediment  concentration  profiles. 

e.  Entrainment  only  through  the  lateral  boundaries  of 
the  jet. 


id 
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Entrainment  Function  and  Similarity  Profiles 


99.  The  quantitative  expression  of  the  entrainment  processes  was 
a  positive  step  toward  the  solution  of  the  turbulent  jet  problem. 


56 


At'.'ZvZ'v'Cv 


Experimentally,.  Ellison  and  Turner  (1959)  found  that  the  entrainment,  E, 
can  be  defined  as  a  function  of  gross  properties  of  the  jet. 


E  =  eu^(Xj)  (25) 

where  u^  =  the  axial  longitudinal  velocity  and  e  =  a  numerical  coeffi¬ 
cient. 

100.  Another  property  of  the  jets  that  simplifies  matters  and 
helps  the  solution  is  similarity.  According  to  this  property,  the 
velocity  profile  remains  similar  to  itself  along  the  various  cross 
sections  of  the  jet.  For  cases  where  bottom  friction  is  important,  the 
similarity  assumptions  must  be  used  with  reservation.  Since  the  exact 
similarity  form  is  not  known,  there  are  a  variety  of  similarity  function 
one  can  choone  from.  These  can  be  either  of  pure  exponential  form  (Fox 
1970)  or  of  pure  polynomial  form  (Stolzenbach  and  Harleman  1971), 

101.  For  this  study  the  similarity  expression,  G(s),  was  chose  to 
lie  between  the  two  aforementioned  cases  (see  Figure  19),  that  is, 

G(s)  =  (1  -  s^)  exp  (-s^)  (26) 

where  s  =  a  nondimensional  parameter. 

102.  The  similarity  assumption  can  also  be  extended  to  the  sedi¬ 
ment  concentration  profiles.  Following  the  suggestion  of  Stolzenbach  and 
Harleman  (1971),  this  similarity  function,  R(s),  can  be  taken  as 

R(s)  =  G(s)‘*  =  (1  -  s^)^exp  (-5ss^)  (27) 

Nondimensional  Form  of  the  Governing  Equations 

103.  Under  all  the  simplifications  and  assumptions  previously 
mentioned,  the  governing  equations  can  be  drastically  reduced  into 
simpler  forms.  Assuming,  furthermore,  that  the  lateral  velocity,  u^,  is 
much  smaller  than  the  longitudinal,  u, ,  and  also  that  the  velocities  are 


vanishing  at  the  side  boundaries  of  the  jet,  the  equations  can  be  treated 
as  follows. 

Continuity  equation 

104.  Under  steady-state  conditions  Equation  10  can  be  written  as 


3(hu  )  9(hu-) 

- —  +  - —  =  0 


dx 


1 


3x, 


(28) 


where  h  is  an  one  dimensional  function  of  the  coordinate. 

105.  Integrating  Equation  28  along  the  width  of  the  jet  yields 


b(xp 


h(x 


l)  J 


X2)dx2 


-b(x,) 


=  2Eh(Xj) 


(29) 


where  b(Xj)  =  half  width  of  the  jet.  A  detailed  description  of  the 
variables  involved  is  given  in  Figure  20. 

106.  According  to  the  similarity  assumption,  the  velocity  is  given 
as 


u^(Xj,X2)  =  u^(Xj)G(s) 


(30) 


where  s  = 


b(Xj) ■ 


107 .  Substitution  of  the  velocity  Equation  30  and  the  entrain¬ 
ment  Equation  25  into  Equation  29  gives 


d 

dx. 


+  1 


h(Xi)Uc(Xi)b(Xi 


)  J  (1  -  s^) 


exp  (-s^)ds 


=  2eh(x^)u^(x^) 


(31) 
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Using  the  Gaussian  quadrature  eight-point  numerical  integration  (see 
Appendix  A) ,  the  evaluation  of  the  integral  was  found  to  be 


(1-s^)  exp  (-s^)  ds  =  1.115 


(32) 


Thus  Equation  31  can  be  written  as 


4 —  (hu  b)  =  Ahu 
dx,  c  c 


(33) 


where  A  =  1.794e.  For  convenience  and  generalization  the  variables  are 
normalized  as 


r  = 


B  = 


^o’ 


H  = 


^o’ 


u  = 


(34) 


where  subscript  0  denotes  the  value  of  the  variables  at  the  outlet  (x^ 
=  0) .  It  can  be  easily  shown  that 


=  AHU 


(35) 


with  H  =  B  =  U  =  1  at  r  =  0. 


(36) 


Equation  35  is  the  nondimensional  continuity  equation  that  constitutes 
part  of  our  mathematical  model. 


Momeptum  equation 

108.  Under  the  same  assinnptions  used  for  the  continuity  equation, 
for  steady-state  conditions  the  expression  for  the  momentum  balance  along 
the  Xj-axis,  Equation  20a  can  be  written  as 


dx 
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b(x^)  2 

h(xp  J  u^  (Xj,X2)dx2 


-b(Xj) 


^  '7 

-b(Xj)-' 


b(Xj)  2 


(37) 


where  f  =  g/C^  =  Darcy-Weisbach' s  friction  coefficient.  Utilizing 
Equation  30,  Equation  37  transforms  into 

+  1 
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h(Xj)u^  (xj)b(x 
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(1  -  s^)^  exp  (-2s^)  ds 


=  -fu^  (Xj) 


+1 
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-1 


(1  -s^)^  exp  (-2s^)  ds 


(38) 


or 


^  (hu  ^b)  =  -fu  ^b 

dXj  c  c 


(39) 


Introducing  one  additional  normalized  variable  for  the  friction  (F)  as 


F  = 


(40) 


and  referring  to  the  normalized  variables  of  Equation  34,  Equation  39 
becomes 


dr 


(HU^B)  =  -FU^B 


(41) 
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subjected  to  the  same  boundary  conditions  given  by  Equation  36. 
Equation  41  is  the  nondimens ional  momentum  equation  which,  along  with 
Equation  35,  constitutes  the  hydrodynamic  part  of  our  mathematical 
model . 

Mass  conservation  equation 

109.  Following  the  same  previous  assumptions  and  neglecting  the 
diffusion  terms  as  being  small  in  comparison  to  the  advective  terms, 
integrating  Equation  24  vertically  over  depth  first  and  then  along  the 
width  of  the  jet  results  in 
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cdx^  =  2Ehc 
2  a 


•/  ■  r 


2  dx2 


where  c  =  sediment  concentration,  c  =  sediment  concentration  in  the 

’  a 

receiving  waters,  w^  =  settling  velocity  of  the  sediment  particles,  and 
u^^  =  critical  velocity  under  which  deposition  occurs.  The  right-hand 
side  of  Equation  42  represents  the  integrated  value  of  S  in  Equation  24. 
More  specifically,  the  first  term  refers  to  the  sediment  gains  due 
to  lateral  entrainment,  while  the  second  term  is  the  sediment  losses  due 
to  deposition.  In  this  model,  no  erosion  processes  are  taken  into 
account.  The  expression  for  the  deposition  function  was  first  introduced 
by  Krone  (1962). 

110.  According  to  the  similarity  assumption,  the  sediment  con¬ 
centration  profiles  are  expressed  as 


c(Xj,X2)  =  c^(Xj)R(s) 


where  c^  =  the  center-line  sediment  concentration.  R(s)  is  given  by 

Equation  27.  Considering  the  receiving  waters  as  sediment-free  (c  = 

a 

0)  and  the  setting  velocity  as  constant  (w^  =  constant),  and  substituting 
Equation  43  for  the  concentration  in  Equation  42,  one  arrives  at 

^[h-^bc^l(|)]  =  -»„bc^l(y  l(|)  (44) 


where 


:(m)  =  J  (1  -  s^) 


™  exp  (-ms^)  ds 


Utilizing  the  Gaussian  quadrature  numerical  integration  of  Appendix  A, 
the  values  of  the  integrals,  I(m),  are  found  to  be 


Ij  =  l(|j  =  0.948;  I2  =  1(1)  =  1.397; 
l3=l(|)=  0.760 


(45a) 


111.  Introducing  the  normalized ^form  of  new  variables  in  addition 
to  those  given  by  Equation  34, 


c  u  b-w- 

c  =  ;  u  w  = 

'0  "  “0  ^“o 


and  substituting  them  into  Equation  44  yields 


^(I,HUBC]  =  -I„WBC  +  I..W  - ;r  C 

dr  1  z  3  y  2 

cr 


Equation  47  can  be  written  as 
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subjected  to  the  boundary  condition 


C  =  1  at  r  =  0 


iLt  t'f  iA  IL.ri 


Equation  48  is  the  nondimensional  fom  of  the  sediment  conservation 
equation  which  along  with  the  hydrodynamic  equations  completes  our 
mathematical  model. 


Analytical  Solutions 


112.  Equations  35,  36,  41,  48,  and  49  formulate  the  mathe¬ 
matical  model  of  a  plane,  sediment-laden  jet  issuing  into  a  shallow 
quiescent  bay.  In  the  following,  solutions  will  be  given  for  the  cases 


of  constant  depth  (H  =  1)  and  of  linearly  varying  depth  (H  =  1  +  arb^/h^, 


a  =  slope).  The  solution  will  be  given  in  terms  of  the  following 
dependent  variables:  velocity  u^Cx^.x^),  width  b(x^),  and  sediment 


concentration  cCx^x^). 


Solution  for  constant  depth  (H  =  1) 

113.  For  constant  depth,  H  =  1,  Equations  35  and  41  become 


^(UB)  =  AU 


—  (U^B)  =  -FU^B 


From  Equation  51  it  is  easily  derived  that 


U  B  =  exp  (-Fr) 


Combining  Equations  52  and  50  thus  eliminating  the  variable  B,  we 


^  +  FU  =  -A  exp  (Fr) 


Equation  53  is  a  Bernoulli-type  ordinary  differential  equation. 

Following  the  typical  solution  procedure  of  Appendix  B,  and  utilizing  the 
boundary  conditions  (Equation  36),  the  nondimens ional  center-line  jet 
velocity  is  given  as 

U  =  ^-  exp  (Fr)  +  exP  (2Fr)J  (54) 

Then  the  solution  for  the  nondimensionalized  jet  width,  B,  can  be 
obtained  directly  from  Equations  52  and  54  as 

B  =  -|^ 

In  dimensional  form  the  jet  velocity,  Uj(r,s),  and  the  jet  width,  b(r), 
can  be  written  as 

Uj(r,s)  =  Uq(1  -  s^)  exp  (-s^)j^-  exp  (Fr) 

-1— ^ 

+  (^1  +  exp  (2Fr)J  (56) 

and 

b(r)  =  b^^-  p  +  (l  +  p) 


114.  For  constant  depth  the  mass  conservation  relation,  Equation 
48,  is  also  reduced  to 

dC  I„  J/  I,  U^\  I,  d(UB) 

—  =  -^(UB)"  -  —  +  - r  WB  -  —  -  dr  (58) 

^  ^1  L\  S  h  J 


dC  I  w  I  WU  d(UB) 

—  - - dr  + - 2  - UB" 

C  I,  U  I,  U 

1  1  cr 
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By  direct  integration  it  results  in 


In  C 


I2  f  1  I3  w  /■ 

=  --  w/  -dr  +  -^—  / 

'10"  ^cr  { 


Udr  -  In  (UB) 

1  “cr  0  (60) 


or 


C  =  exp 


I2  f  1  h  ^  f 

- W  ;  -  dr  +  -^  - -  /  Udr  -  Jin  (UB) 

I,  U  I  U 

1  0  ^1  ^cr  0 


(61) 


The  integrals  within  the  exponential  cannot  be  evaluated  analytically. 
Thus,  for  the  application  of  Equation  61,  a  numerical  integration 
technique  is  required. 

115.  An  essential  simplification  in  the  solution  can  be  achieved 
by  making  an  additional  assumption  of  no  entrainment  conditions  (A  =  0) 
in  Equation  54.  Then  the  integrals  can  be  easily  evaluated  as 


and 


/ 

/ 


exp  (Fr)  dr  = 


exp  (-Fr)  dr  = 


I  exp  (Fr) 


1 

-f  exp 


(-Fr) 


(62) 


(63) 


Therefore  Equation  60  can  be  rewritten  as 


I2  w 


I3  W 


In  C  = - exp  (Fr)  - - —  exp  (-Fr)  +  Y 

T  T?  r  11  fci-. 


(64) 


I,  U  “F 
1  cr 


Since  An(UB)  -  £n[exp(-Fr)exp(Fr) J  -  £n(l)  =  0,  from  the  boundary  con¬ 
dition,  Equation  49,  the  constant  of  integration,  Y,  is  found  to  be 


Thus,  finally,  the  solution  for  the  normalized  center-line  sediment 
concentration  reads 


I  W  I  W  1 

C  =  exp  ~  [1  ■  exp  (Fr)  ]  +  — - ^  [1  -  exp  (-Fr)]  ^  (66) 


I-  F  U  ‘ 
1  cr 


where  ~  ~  0-801*  1°  terms  of  the  independent 

variables,  r  and  s,  the  sediment  concentration  within  the  jet  is  given 


c(r,s)  =  Cq(1  -  s  exp  ( 


/12\  (^2'^ 


=  [1  -  exp  (Fr)] 


I-  W  1  ^ 

+ - 2  [1  -  exp  (-Fr)  1  ((67) 

I,  F  U  I 

1  cr  / 


Solution  for  linearly  varying  depth  (H  =  1  arb^/hQ) 

116.  For  linearly  varying  depth  Equations  35  and  41  are 
written  as 


d(UB)  b 

H  — : -  +  a  r^UB  =  AHU 

dr  hg 


(68) 


and 


d(U‘B)  /  b  \ 

H  — =  -K  ^7^ 


(69) 


Solving  Equation  69  for  U  B  gives 


where  D  =  f/a.  The  combinatiOQ  of  Equations  70  and  68  yields 


Hj-  +  a^  =  AHU 


+  a^(D  -  1)h"^(HU)  =  -AH®‘\hU)^ 


Equation  72  is  also  a  Bernoulli-type  ordinary  differential  equation 
with  respect  to  the  monomial  HU.  Thus,  along  with  the  boundary  condi¬ 
tions  of  Equation  36,  the  solution  for  HU  (Appendix  B)  is 

HU  =  ^  54-d  +  S  (73) 


Consequently,  the  normalized  center-line  velocity  is 


u  =  caArV"  ^  5-1^  (H^-o -1)  .  ^ 


The  normalized  jet  width  can  now  be  derived  from  Equations  74  and  70 


®  abg  2  -  D  ^  2A 


Following  all  these,  the  velocity  and  the  width  of  the  jet  in  terms  of 
the  variables  r  and  s  are  given  as 


Uj(r,s)  =  Uq(2A)-V°  ^ 


h„  1 


2  -  D 


(H^"°  -1)  + 


-T' 

2A 


.  (1  -  s^)  exp  (  -s^)  (76) 
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and 


b(r)  =  bQ2AH°‘ 


ab. 


2  -  D 


(H 


2-D 


•1)  ^2A 


(77) 


117.  For  linearly  varying  depth,  the  sediment  conservation 
relation,  Equation  48,  after  direct  integration,  reads 


or 


In  C 


r 


r 

U 

-  dr  -  £n  (HUB)  (78) 
H 

0 


I 


I 


C  =  exp 


+ 


r 


/ 


U 

-dr  -  £n(HUB) 
H 


0 


(79) 


Again,  the  integrals  within  the  exponential  must  be  evaluated  numeri¬ 
cally.  In  terms  of  the  variables  r  and  s,  the  sediment  concentration  is 
written  as 


c(r,s)  =  CqC(1  -  s^)"*  exp  y-  |  s^j  (80) 

where  C  is  given  by  Equation  79. 

Approximation  of  Sediment  Deposition  in  a  Bay 

118.  The  sediments  carried  by  the  river  waters  are  the  main 
material  for  the  river  delta  development.  The  form  of  the  delta, 
however,  depends  on  other  parameters  such  as  the  climate,  the  hydro- 
dynamic  field,  the  area  topography,  and  the  sediment  characteristics 
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themselves.  Neglecting  erosion  and  side  entrainment  of  the  sediment, 
under  steady-state  conditions,  material  deposition  is  linear  in  time  and 
can  be  computed  as 


where  c^  =  deposited  sediment  and  t  =  time.  The  values  for  the  velocity, 
Uj,  and  the  concentration,  c,  are  taken  respectively  from  Equations 
56  and  67  or  76  and  80. 

119.  In  the  case  that  the  variable  c  is  given  in  units  of  mass  per 
volume,  the  deposited  sediment,  c^,  is  computed  in  mass  per  unit  area. 
Therefore,  assuming  no  consolidation  processes,  the  thickness  of  the 
deposited  sediment  layer,  d,  can  be  directly  computed  from  Equation  81 
as 


d 


(82) 


where  =  density  of  the  sediment. 

120.  For  the  numerical  simulation  of  river  delta  evolution,  the 
overall  phenomenon  is  taken  as  steady  for  certain  time  intervals.  Then, 
due  to  the  variation  in  the  bottom  topography  from  deposited  sediment 
buildup,  proper  adjustments  must  be  made  to  various  parameters  and  the 
solution  repeated  under  the  new  conditions.  In  such  a  manner  the  mathe¬ 
matical  model,  although  being  of  a  steady-state  nature  itself,  predicts 
the  delta  development  over  time. 
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PART  VI:  PREDICTION  OF  THE  ATCHAFALAYA  RIVER  DELTA  GROWTH 


Basic  Data  for  Analytical  Prediction 

121.  The  Atchafalaya  River-Bay  area  is  a  complicated  dynamic 
physical  system.  For  the  purpose  of  this  study,  the  values  of  various 
factors  controlling  this  system  must  be  known.  The  most  essential  of 
these  factors  are  the  geometry  and  topography  of  the  area;  the  water 
discharge,  and  sediment  load  of  the  river;  and,  to  a  lesser  degree,  the 
hydrology  and  climate  of  the  surrounding  environment.  Field  data  per¬ 
taining  to  our  analytical  analysis  of  delta  development  were  gathered 
from  various  sources  and  are  outlined  in  the  following  sections. 
Physical  dimensions  of  river  outlets 

122.  Two  river  outlets,  the  Lower  Atchafalaya  River  Outlet  and 
Wax  Lake  Outlet,  are  forming  in  Atchafalaya  Bay.  Each  outlet  has  its 
own  characteristics  and  these  consequently  influence  the  form  of  the 
developing  delta. 

123.  Using  data  from  a  1974-76  hydrographic  analysis  by  the 
Waterways  Experiment  Station  (1981),  the  cross  section  of  the  Lower 
Atchafalaya  River,  at  river  mile  135.8  from  Simmesport,  La.,  is  plotted 
in  Figure  21a.  The  river  outlet  is  approximately  3,000  ft  wide  and 
reaches  a  depth  of  35  ft.  From  the  same  data  source,  the  cross 
section  of  Wax  Lake  Outlet,  at  river  mile  122.3  from  Simmesport,  is 
plotted  in  Figure  21b;  the  outlet  is  about  1,000  ft  wide  and  reaches  a 
depth  of  56  ft. 

124.  The  Atchafalaya  River  Outlet  is  under  continuous  dredging 
so  that  a  navigation  channel  can  be  maintained  throughout  the  bay. 

This  navigation  channel  is  developing  as  the  future  main  course  of  the 
river  through  the  emerging  delta. 

Initial  bathymetry  conditions  of  the  receiving  bay 

125.  Bathmetric  maps  of  the  Atchafalaya  Bay  for  the  years  1972 
and  1977  (Adams  and  Baumann  1980)  give  a  clear  picture  of  bay 
bathymetry  (Figures  8  and  13).  With  a  depth  ranging  from  4  to  8  ft 
(1.2  to  2.4  m) ,  the  bay  can  be  classified  as  shallow  and  well-mixed. 


thus  bottom  friction  may  play  a  predominant  role.  Other  indications 
from  these  bathymetric  maps  are  that  the  newly  deposited  sediment  is 
close  to  the  river  mouth,  and  that  the  bathymetry  of  the  remaining  bay 
is  less  affected. 

River  discharge  and  suspended  sediment  load 

126.  The  water  discharge  of  the  Atchafalaya  River  at  Simmesport 
(Figure  11)  and  at  the  Lower  Atchafalaya  River  Outlet  (Figure  12)  is 
varied.  It  ranges  from  80,000  cfs  to  600,000  cfs  (2,000  cms  to 
17,000  cms). 

127 .  The  current  velocity  at  the  river  mouth  ranges  from  3  fps 
(1  mps)  to  0.2  fps  on  a  diurnal  basis  due  to  the  tidal  action  (Figure 
22).  However,  the  current  always  has  a  southwest  direction  at  the 
Lower  Atchafalaya  River  Outlet  (WES  Preliminary  Field  Data  Report 
1982). 

128.  The  character  of  the  sediments  entering  the  bay  via  riverine 
flows  is  controlled  by  the  sedimentation  processes  within  the  river-bay 
system.  According  to  Roberts,  Adams,  and  Cunningham  (1980),  the  sediment 
load  is  composed  mainly  of  silt  and  clay  (>  75  percent)  and  a  small  part 
of  sand  (<  25  percent)  (Table  2).  Thus  our  study  can  be  simplified  by 
dealing  with  the  fine  sediment  materials  which  are  transported  in  sus¬ 
pension  only.  Whether  to  treat  the  sediment  as  cohesive  or  noncohesive 
is  a  difficult  question.  Since  much  of  the  subaerial  sediment  is  sandy, 
it  is  reasonable  to  use  a  single  equation  to  predict  sedimentation 
processes  without  considering  various  particle  sizes. 

129.  Letter  (1982),  based  on  observations  of  the  sediment  con¬ 
centrations  in  Simmesport,  derived  a  regression  equation  relating 
water  and  sediment  discharge  obtaining 

Q  =  0.0728  Q  (83) 

s  w 

where  is  the  suspended  sediment  load  in  1000  tons  (2000  lb)  per 
day  and  is  the  water  discharge  in  1000  cfs.  The  value  of  the 
exponent  is  close  to  unity,  permitting  the  assumption  of  a  linear 
relationship  between  water  and  sediment  discharge. 
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130.  As  shown  in  Figure  12,  the  variations  of  water  and  sedi¬ 
ment  discharge  at  the  Lower  Atchafalaya  River  Outlet  follow  the  same 
pattern  through  time.  Thus,  for  all  practical  purposes,  the  mean 
sediment  concentration  can  be  reasonably  assessed.  A  mean  water 
discharge  of  300,000  cfs  (8,500  cms)  and  a  mean  suspended  sediment 
load  of  200,000  metric  tons  per  day  are  shown  by  Figure  12  for  the 
period  1973-75.  From  these  values,  an  average  sediment  concentration 
of  0.27  kg/m  (270  ppm)  is  derived. 

131.  For  a  water  discharge  of  300,000  cfs.  Equation  83  gives 

3 

a  mean  sediment  concentration  of  275,000  tons/day  (0.34  kg/m  ),  which 

3 

is  close  to  the  value  of  0.27  kg/m  .  In  the  following  computations, 
the  mean  sediment  concentration  (Cq),  varying  from  100  ppm  to  600 
ppm,  is  taken  into  consideration  for  our  study. 

Sediment  settling  velocity  and  deposition 

132.  The  sediment  deposition  rate,  which  depends  on  the  shear 
stress  of  the  flow  (Krone  1962),  is  assumed  to  be  proportional  to  the 
velocity  of  the  sediment  particles.  In  general,  this  velocity  depends 
on  the  shape,  size,  and  weight  of  the  particle,  as  well  as  on  hydro- 
dynamic  conditions.  To  define  the  settling  velocity  of  cohesive 
sediments  in  a  real  situation  is  a  difficult  task,  because  these 
particles  have  the  tendency  to  adhere  to  each  other  and  form  large 
aggregated  floes.  Laboratory  measurements  made  by  the  Waterways 
Experiment  Station  on  Atchafalaya  River  sediments  (Figure  23)  have 
shown  a  settling  velocity  from  0.01  mm/sec  to  <  1.0  oun/sec  (WES 
Preliminary  Field  Data  Report  1982). 

133.  The  freshly  deposited  sediments  are  initially  in  a  very 
loose  state.  However,  as  new  sediment  layers  are  superimposed,  the 
low  density  large  aggregates  are  crushed  down  to  smaller  floes,  and 
denser  and  stronger  bed  layers  are  formed.  As  a  result  of  this 
consolidation  effect,  the  bulk  density  of  the  bed  (p^)  may  be  of  high 
value  (Krone  1978).  Wells  and  Roberts  (1980)  used  a  density  of  375 

3 

kg/m  to  estimate  the  sediment  transported  in  the  Atchafalaya  mud 
stream.  In  our  analysis,  a  value  of  equal  to  400  kg/m^  is  used. 


76 


STA. 'H;  10/29/01;  S=0. 3ppt;  Co=754pp 


e.  1000 


134.  The  active  development  of  the  river  deltas  in  Atchafalaya 
Bay  indicates  that  the  deposition  rate  in  the  area  is  much  faster  than 
the  rate  of  erosion.  Normally,  the  combined  effect  of  both  deposition 
and  erosion  should  be  considered.  In  our  analytical  model  it  is 
assumed  that  the  sediment  is  subjected  only  to  transportation  and 
deposition,  and  that  each  particle  that  hits  the  bottom  adheres  to  it. 
Supressing  the  resuspension  effects  is  acceptable,  since  the  sediment 
deposition  is  estimated  based  on  average  values  of  the  water  discharge 
and  the  sediment  concentration. 

Bottom  resistance  and  lateral  entrainment 

135.  The  predominance  of  silt,  clay,  and  fine  sand  in  the 
bottom  materials  of  the  bay  causes  the  bed  surface  to  be  smooth,  with 
relatively  small  friction  resistance.  In  the  absence  of  field  data, 
the  Darcy-Weisbach  coefficient  of  friction  (f)  is  assumed  to  vary 
from  0.001  to  0.006  in  this  study.  It  is  further  assumed  that  the 
friction  coefficient  is  time- independent . 

136.  The  riverine  waters  entering  Atchafalaya  Bay  are  well 
mixed  with  the  receiving  waters.  Because  of  the  shallowness  of  the 
bay,  there  is  no  density  stratification,  and  the  only  entrainment  of 
the  receiving  bay  waters  to  the  sediment-laden  river  waters  is  through 
the  sides  of  the  jet.  The  lateral  entrainment  coefficient  (e)  has 
proven  to  be  a  function  of  Richardson's  number  (Ellison  and  Turner 
1959).  Engelund  (1976)  used  a  value  of  0.075  in  the  case  of  a  very 
small  Richardson  number.  In  our  analysis,  a  range  of  e  varying  from 
0.0375  to  0.300  is  considered. 


Procedures  for  Closed-Form  Analytical  Solutions 


Idealization  of  the  Atchafalaya  River-Bay  system 

137.  The  physical  features  of  the  Atchafalaya  River-Bay  system 
are  approximated  by  simple  geometry  for  the  domain  of  the  closed-form 
analytical  solution.  Using  the  data  of  the  previous  section,  the 
nominal  values  of  dependent  variables  (h,  b,  u,  and  c)  can  be  deduced 
for  our  study. 


138.  The  cross  sections  of  the  Lower  Atchafalaya  River  and  Wax 

2 

Lake  Outlet  (Figures  21a  and  Zlb)  are  found  to  be  equal  to  66,000  ft 
2  2  2 

(6,200  m  )  and  29,000  ft  (2,700  m  ),  respectively.  Considering  a 
mean  annual  water  discharge  of  300,000  cfs  (8,500  cms)  at  Simmesport 
and  the  70-30  percentage  split  of  discharge  (Letter  1982) ,  the  flows 
in  the  Lower  Atchafalaya  River  Outlet  and  the  Wax  Lake  Outlet  are 
210,000  cfs  and  90,000  cfs,  respectively.  Both  flows  exhibit  a 
velocity  of  3.10  cfs  (^  1  mps)(WES  1981). 

139.  Basically,  the  river  mouth  is  the  point  at  which  fresh 
water  leaves  the  confined  channel  and  mixes  with  ambient  water.  In 
the  Atchafalaya  River-Bay  system,  the  river  outlet  discharges  a  part 
of  its  flow  into  the  bay  through  a  navigation  channel.  The  bottom  of 
the  navigation  channel  is  much  deeper  than  the  bottom  of  the  ambient 
bay.  In  this  study,  those  flow  and  sediment  discharges  through  the 
navigation  channel  are  assumed  not  to  have  significant  influence  on 
jet  characteristics,  and  are  assumed  to  deliver  directly  to  the  Gulf 
of  Mexico. 

140.  The  Atchafalaya  Bay  is  relatively  flat,  with  a  uniform 
depth  of  6.0  ft  (~2  m)  (Figure  13).  This  depth,  measured  below  mean 
sea  level,  is  referenced  to  the  bathymetry  conditions  of  the  year  1977 
(Adams  and  Baumann  1980).  In  a  shallow  bay,  the  formation  of  delta 
lobes  depends  on  the  bathymetry  of  the  receiving  bay  rather  than  the 
geometry  of  the  river  itself  (Wells,  Chinburg,  and  Coleman  1984).  This 
fact  leads  us  to  assume  that  the  river  outlet  can  be  approximated  by  a 
rectangular  cross  section.  Using  the  bathymetry  map  of  Adams  and  Baumann 
(1980),  the  inferred  width  is  about  4,000  ft  (1,200  m)  for  the  Lower 
Atchafalaya  River  Outlet,  and  about  3,300  ft  (1,000  m)  for  the  Wax 

Lake  Cutlet.  The  nominal  water  depth  at  both  outlets,  h^  =  6.5  ft 
(2  m) ,  is  used  in  this  study. 

Values  of  various  parameters  used  in  the  analytical  study 

141.  As  previously  mentioned,  the  rate  of  sediment  deposition 
depends  on  the  shear  stress  of  the  flow  (Krone  1962).  In  this  study, 
a  linear  relationship  between  the  shear  stress  and  mean  square 
velocity  is  assumed  (Equation  81).  The  center-line  velocity  at  the 


river  outlet  is  used  as  a  critical  velocity  (u  =  u„)  under  which 
deposition  occurs.  Equation  81  implies  that  at  the  middle  point  of 
the  river  outlet  no  deposition  occurs.  Table  4  summarizes  the  range 
of  the  values  for  various  parameters  to  be  used  in  our  analytical 
study. 


Computer  graphics  and  numerical 
integration  used  for  analytical  solutions 

142.  The  SAS/GRAPH  (1981),  a  computer  graphic  system,  is  used 
to  display  the  values  of  dependent  variables  (u,  b,  c,  d)  for  various 
cases.  The  GPLOT  PROCEDURE  graphs  one  variable  against  another, 
producing  a  two-dimensional  plane.  The  G3D  PROCEDURE  plots  the  value 
of  three  variables  and  produces  a  three-dimensional  surface.  The 
variables  to  be  plotted  are  specified  in  a  PLOT  statement.  Both  the 
GPLOT  and  G3D  PROCEDURES  can  automatically  scale  the  axes,  or  the  user 
can  specify  the  scale.  Use  of  computer  graphics  provides  flexibility 
in  displaying  information  meaningfully.  Examples  for  use  of  computer 
graphic  programs  are  listed  in  Appendix  C. 

143.  A  numerical  integration  technique  is  needed  for  t.  le  compu- 
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tation  of  integrals,  dr  and  Jq  dr,  appearing  in  Equation  79. 

The  orthogonal  collocation  method  is  used  for  numerical  integration 

(Villadsen  and  Michelson  1978).  The  method  consists  of  expanding  the 

normalized  center-line  velocity  U,  a  dependent  variable,  in  terms  of  a 
t  h 

n  order  Jacobi  polynomial.  The  n  roots  of  this  orthogonal  poly¬ 
nomial  are  chosen  as  the  n  collocation  points  (Kuu  and  Polack  1982). 
The  integration  of  the  profile  of  the  dependent  variable  is  approxi¬ 
mated  by  the  Radau  quadrature  formula  (Villadsen  and  Michelsen  1978). 
The  Radau  quadrature  weights  at  the  n  collocation  points  (the 
abscissas)  are  determined  (see  Appendix  D) .  The  integral  is  approxi¬ 
mated  by  the  summation  of  the  product  of  the  weight  and  the  value  of 
the  function  at  the  collocation  points.  Quadrature  integration  is 
explained  in  detail  in  Appendix  D. 


Base  Results  of  Analytical  Solutions 


144.  This  section  presents  a  basic  computation  using  nominal 
values  of  variables  to  obtain  the  closed-form  analytical  solution. 

This  base  result  will  serve  as  a  guideline  for  the  comparison  of  the 
results  of  various  cases  and  for  the  sensitivity  analysis  of  various 
parameters . 

145.  The  following  nominal  values  are  used  in  the  basic  computa¬ 
tions  to  obtain  the  base  results: 
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Results  of  constant  depth 

without  entrainment  (Case  1:  E  =  0,  H  =  1) 

146.  The  normalized  (dimensionless)  jet  center-line  velocity  (U) 
and  jet  width  (B)  are  obtained  by  substituting  A  =  0  into  Equations 

81 


54  and  55.  The  center-line  concentration  (C)  is  given  in  Equation 
66.  The  equations  for  each  of  these  variables  are 


U  =  exp  (-Fr) 


B  =  exp  (Fr) 


i  I  W  I  W  1 

C  =  exp/  - - (1  -  exp  (Fr)]  +  — - :r  [1  -  exp  (-  Fr)] 


I,  F  U 
1  cr 


147.  The  dimensional  form  of  jet  velocity  (u) ,  jet  width  (b) , 
and  center-line  concentration  (c)  are  derived  straightforwardly  as: 


u(r,  s)  =  u^G(s)  =  UqG(s)U 


=  u-(l  -  s^)exp  (-s^)U 


b(r)  =  b^B 


c(r,  s)  =  c^R(s)  =  CqR(s)C 


=  Cq(1  -  s^)"*  exp  (-^s^)C 


The  deposited  sediment  (c^)  and  its  thickness  (d)  are  computed  from 
Equations  81  and  82,  respectively: 

Cd(r,  s)  =  Woc(l  -  u^/u^/)t  (9 


d(r,  s)  -  c^/Pg 


The  results  of  the  calculations  of  u,  c,  and  d  are  displayed  three- 
dimensionally  in  Figures  24,  25,  and  26,  respectively  (and  Plate  1), 
with  the  aid  of  computer  graphics  (see  Appendix  C).  Each  of  the 
variables  is  plotted  versus  the  normalized  longitudinal  (r)  and 
lateral  (s)  coordinates. 
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Figure  25.  Case  1:  Base  results  for  sediment  concentration  c 
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Results  of  linearly  varying  depth 
without  entrainment  (Case  2;  E  =  0,  H  ^  1) 

148.  For  linearly  varying  depth  (H  =  1  +  arb^/h^)  and  no 

entrainment  (A  =  0),  the  nondimens ional  continuity  equation.  Equation 

35,  becomes 

—  (HUB)  =  0 
dr 


Using  the  boundary  conditions  stated  in  Equation  36,  Equation  92 
implies  that 


HUB  =  constant  =  1 
or 

UB  =  H'^ 


(93) 

(94) 


Similarly,  the  nondimensional  momentum  equation,  Equation  41,  can 
be  written  as 


d  2  ^9 

H  —  (U^B)  +  (a  —  +  F)U^B  =  0 

dr  h^ 


(95) 


Solving  for  U  B,  we  obtain 


U^B  =  H-‘'  "  “> 


Combining  Equations  94  and  96,  we  get 


U  =  H 


B  =  H 


■1  +  D 


and  finally 


(96) 


(97) 

(98) 
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HUB  =  (H)  (H  )  (H  *  )  =  1  (99] 

149.  The  nondimensional  form  of  sediment  conservation  is  given 
by  Equation  48.  Substituting  Equations  97,  98,  and  99  into 
Equation  48,  we  obtain 


+  —  - 5 

“cr  / 


C  =  exp 


f  Vl  ^  “cr  /  . 


where 


Y  =  -  ^ 

F  Vl, 


C  =  exp 


■I2  w 

.'1  " 


I  W 

(1  -H°) 


(1  - 


I,  F  U  ‘ 
1  cr 


150.  The  dimensional  form  of  jet  velocity  (u) ,  width  (b) , 
concentration  (c),  the  deposited  sediment  (c^),  and  the  sediment 
thickness  (d)  for  Case  2  are  computed  using  exactly  the  same  procedure 
as  Case  1  with  Equations  87,  88,  89,  90,  and  91.  The  three- 
dimensional  graphs  for  u,  c,  and  d  for  this  case  are  plotted  in 
Plate  2. 

Results  of  constant  depth 
with  entrainment  (Case  3:  E  ^  0,  H  =  1 


151.  As  shown  in  PART  V  and  presented  again  in  this  section, 
Equations  54,  55,  and  61  are  the  expressions  for  the  normalized 
center-line  velocity  (U) ,  jet  width  (B) ,  and  center-line  concentration 
(C),  respectively: 


87 


(104) 


"=[-r 


exp  (Fr)  +  (1  +  exp  (2Fr) 


B  =  -  p  exp  (Fr) 


(105) 


C  =  exp 


r  ^2  r  ^  h  ^  r 

-  —  W  /  —  dr  +  —  - 5^  / 

L  1  0  1  cr  0 


Udr  -  £n(UB) 


(106) 


152.  The  integrals  appearing  in  Equation  106  have  to  be  eval¬ 
uated  numerically.  In  this  case,  the  lower  end  point  (U  =  1  at  r  =  0) 
of  both  integrals  is  known,  thus  the  Radau  quadrature  formula 
(Appendix  D)  is  used. 

153.  Again,  the  dimensional  form  of  jet  velocity  (u) ,  width 


(b) ,  concentration  (c) ,  the  deposited  sediment  (c^)  and  its  thickness 


(d)  are  computed  straightforwardly  by  using  Equations  87,  88, 
89,  90,  and  91,  respectively.  Plate  3  contains  the  plots  for 
u,  c,  and  d. 


Results  of  linearly  varying  depth 
with  entrainment  (Case  4:  E  ^  0,  H  /  1) 


154.  The  general  solution  of  linearly  varying  depth  with 
entrainment  for  the  normalized  center-line  velocity  (U) ,  jet  width  (B) , 
and  center-line  concentration  (C)  are  given  in  Equations  74,  75, 


and  79  in  PART  V: 


U  =  (2A)-'‘  H-"  ^  ^ 


(107) 


B  .  :-5-  ^  ^ 


(108) 


r 


and 


C  =  exp 


U 

-dr  -  iln(HUB) 
H 


(109) 


155.  Again,  the  lowest  limit  (H  =  1,  U  =  1,  at  r  =  0)  of  both 
integrals  appearing  in  Equation  109  is  known,  and  the  Radau  inte¬ 
gration  technique  (Appendix  D)  is  used.  Similarly,  the  dimensional 
form  of  jet  velocity  (u) ,  jet  width  (b),  concentration  (c),  the 
deposited  sediment  (c^),  and  the  sediment  thickness  (d)  are  computed 
by  using  Equations  87,  88,  89,  90,  and  91,  respectively.  The  three- 
dimensional  graphs  of  u,  c,  and  d  are  displayed  in  Figures  27,  28, 
and  29,  respectively  (and  Plate  4). 


Prediction  of  Delta  Front  Advancement 

156.  The  rate  of  delta  growth  depends  on  the  amount  of  sediment 
supplied  by  the  riverine  waters  and  reworking  by  current  forces  in  the 
receiving  bay  (Coleman  and  Wright  1975).  The  areal  and  mass  extent 

of  deltaic  evolution  is  governed  by  the  relative  roles  of  inertial  and 
frictional  forces.  Thus  sediment  deposition  patterns  are  determined 
by  various  physical  parameters  that  are  formulated  in  the  analytical 
solutions  for  the  four  different  cases. 

Sediment  deposition  patterns  under  quasi-steady  state 

157.  In  this  study,  a  quasi-steady  state  for  sediment  deposi¬ 
tion  is  assumed  (Engelund  1976).  The  deposited  sediment  and  its 
thickness  are  computed  as  a  linear  function  with  time  (Equation  90). 
The  deposition  patterns  of  Case  4  for  the  time  interval  of  1.5  years, 
1.0  year,  and  0.5  year  are  given  in  Plate  5.  From  these  figures,  it 
is  shown  that  the  deposited  sediment  forms  a  saddle-shaped  bottom. 

The  rapid  accumulation  of  suspended  sediment  near  the  river  outlet  and 
the  abrupt  decline  of  sediment  deposition  away  from  the  outlet  are 
observed.  As  the  central  portion  of  the  sediment  accumulates,  it 
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Figure  28.  Case  4:  Base  results  for  sediment  concentration  c 
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causes  the  river  mouth  channel  to  separate  into  arms  known  as  bifur¬ 
cating  channels  (Coleman  1976). 

158.  Therefore,  until  subaerial  land  has  emerged,  it  can  be 
inferred  that  the  riverine  input  has  limited  influence  within  a 
certain  longitudinal  distance;  and  that  as  the  deposition  process  is 
completed  in  this  level  of  development  the  river  mouth  advances  to 
the  end  of  the  subaerial  land  and  begins  the  process  of  bifurcation. 
Conceptualization  of  delta-channel  development 

159.  The  geometry  of  river-mouth  sandbars  is  determined  by 
riverine  flow  conditions.  When  river  outflow  velocities  are  high  and 
water  depths  seaward  of  the  mouth  are  shallow,  the  rapid  rate  of 
effluent  expansion  provides,  initially,  a  broad  radial  sandbar  and 
later  on  develops  a  distributary  network.  Coleman  (1976)  documented 
three  major  types  of  existing  delta  channel  patterns  (Figure  30).  In 
an  environment  having  a  high  subsidence  rate,  low  wave  and  tide 
energy,  low  offshore  slope,  and  a  fine  grain  sediment  load,  the 
development  of  bifurcating  channels  is  typical.  Deltas  developing 
this  distributary  pattern  are  characterized  by  a  large  number  of  river 
mouths  (Figure  30). 

160.  Channel  bifurcations  and  crevasse  discharges  have 
influenced  the  shape  of  subaerial  land  in  Atchafalaya  Bay.  Adams  and 
Baumann  (1980)  identified  five  levels  of  bifurcation  that  have  been 
taking  place  in  the  Lower  Atchafalaya  River  delta  (Figure  31).  The 
branching  channels  discharge  water  and  suspended  sediments  that  form 
subdeltas  along  the  distributary  channels.  The  shoaling  at  the 
branching  river  mouth  is  repeated  and  new  bifurcations  develop  at  the 
new  channel  outlets  to  form  a  complex  branching  pattern  (Adams  ar... 
Baumann  1980) . 

Stepwise  procedure  for  delta  growth  prediction 

161.  A  stepwise  procedure,  formulated  to  estimate  the  areal  and 
volume  extent  of  the  delta  in  both  space  and  time,  is  the  means  by 
which  delta  growth  is  predicted.  The  space  steps  are  selected  at  the 
beginning  of  the  computation  and  guided  by  the  dimensional  plot  of 
sediment  deposition.  The  time-steps  are  not  known  a  priori,  but 


depend  upon  the  analytical  solutions  themselves  and  are  obtained  by 
the  numerical  procedure  (Appendix  E). 

162.  In  order  to  simulate  the  process  of  bifurcation,  it  is 
assumed  that  at  each  level  of  bifurcation  the  branching  channels 
achieve  a  stage  of  development  similar  to  their  parent  channel,  and 
that  the  subdelta  is  produced  by  turbulent  jets  with  the  same  riverine 
flow  conditions.  An  idealized  bifurcation  scheme  is  displayed  in 
Figure  32. 

163.  In  this  study,  our  predictions  are  based  on  mean  sea 
level,  which  is  used  as  the  determining  elevation  for  subaerial  land. 
Plate  6  is  a  two-dimensional  plot  of  sediment  thickness  for  Case  4. 

The  total  sediment  volume  is  computed  by  analytical  integration.  The 
time  required  to  fill  the  known  volume  of  sediment  to  an  average 
thickness  h  is  obtained  by  the  method  of  Bisection  Search  (Scheid 
1968).  The  numerical  procedures  are  explained  in  detail  in  Appendix  E. 
Prediction  for  the  Lower  Atchafalaya  River  Outlet  delta 

164.  River  delta  development  depends  on  the  quantity  of 
sediment  that  is  delivered  to  and  retained  in  the  receiving  bay.  In 
Atchafalaya  Bay,  much  of  the  flow  through  both  the  Lower  Atchafalaya 
River  Outlet  and  Wax  Lake  Outlet  enters  the  Gulf  of  Mexico  via  the 
navigation  channels.  Adams  and  Baumann  (1980)  estimated  that  40 

to  50  percent  of  the  total  volume  of  suspended  sediment  delivered 
to  the  bay  through  the  river  outlets  will  be  retained  in  Atchafalaya 
Bay  and  the  rest  will  be  dispersed  in  peripheral  marshes  and  offshore 
regions. 

165.  For  predictive  purposes,  it  is  assumed  that  the  volume  of 
sediment  is  proportional  to  the  volume  of  discharge.  Based  upon  the 
70-30  percentage  split  of  discharge  for  the  two  outlets  (Letter 
1982),  the  Lower  Atchafalaya  River  delta  is  expected  to  grow  to  a  much 
greater  extent  than  the  Wax  Lake  Outlet  delta. 

166.  The  Lower  Atchafalaya  River  delta  can  be  divided  into  an 
eastern  and  a  western  component  (Adams  and  Baumann  1980),  areas  that 
are  associated  with  East  Pass  and  West  Pass  (Figure  31).  Van  Heerden 
(1980)  conducted  an  extensive  field  study  in  the  eastern  half  of  the 
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Lower  Atchafalaya  River  delta  from  1973  to  1979.  The  total  area  of 

2  2 

the  eastern  half  of  the  subdeltas  was  4.56  mi  (11.67  km  )  with  an 

2  2 

average  growth  rate  of  0.76  mi  /year  (2  km  /year). 

167.  In  contrast,  subdelta  development  on  the  western  side  is 
more  complex.  Adams  and  Baumann  (1980)  indicated  that  if  it  were  not 
for  the  navigation  channel,  God's  Pass  and  Log  Island  Pass  (Figure 
31),  which  developed  from  second-order  bifurcation,  along  with  East 
Pass  would  represent  the  three  major  river  outlets. 

168.  To  simulate  the  growth  of  river  deltas  in  this  physical 
domain,  a  stepwise  procedure  described  in  the  previous  section  is  used 
and  the  simulation  is  done  as  follows: 

a.  In  the  first  order  bifurcation,  the  total  river  discharge 
and  suspended  sediments  are  divided  into  three  equal 
amounts  representing  the  three  major  outlets. 

b.  A  space  step  is  selected  at  the  onset  of  computation;  the 
subdelta  area  (jet  length  and  width)  at  each  outlet  is 
calculated. 

£.  The  time-step  is  searched  by  the  numerical  procedure 
(Appendix  E). 

The  numerical  procedure  is  repeated  in  the  same  manner  for  the  se¬ 
quential  order  of  bifurcation  processes.  The  results  of  each  computa¬ 
tion  are  summarized  in  Table  5.  The  average  growth  rate  is  about  5 
km^/year . 

Prediction  for  the  Wax  Lake  Outlet  delta 

169.  The  bifurcation  process  is  not  as  evident  on  the  Wax  Lake 
Outlet  delta.  Most  of  the  delta  development  has  occurred  west  of  the 
main  channel  (Adams  and  Baumann  1980).  A  few  branching  channels  have 
developed  on  the  western  side  of  the  outlet  and  subdeltas  have  formed 
along  the  subchannels. 

170.  Wells,  Chinburg,  and  Coleman  (1984)  reported  that  during  the 

1980-81  flood  year,  the  Wax  Lake  Outlet  delta  represented  17  percent  of 

total  subaerial  land  in  Atchafalaya  Bay,  and  approximately  10  percent  of 

the  total  if  averaged  over  a  6-year  period  from  1975  to  1981  covering  an 
2  2 

area  of  2  km  to  3  km  .  The  average  growth  rate  ranged  from  0.3  to  0.5 
km  /year. 


171.  To  simulate  the  natural  bifurcation  pattern  in  the  absence 
of  field  observation  is  a  difficult  task.  Thus  the  same  branching 
scheme  presented  earlier  is  used  for  modeling  the  growth  of  the  Wax 
Lake  Outlet  delta.  The  stepwise  numerical  procedure  is  again  used  to 

estimate  the  growth  of  delta  lobes.  The  results  are  listed  in  Table  6. 

2 

The  average  growth  rate  is  about  2.6  km  /year. 

Estimation  of  Atchafalaya  River  delta  growth 

172.  Aerial  photographs  (Wells,  Chinburg,  and  Coleman  1984)  and 
photomosaics  (Adams  and  Baumann  1980)  show  that  the  Atchafalaya  River 
deltas  have  grown  by  developing  parabolic  lobes  of  fine-grain  sediments 
that  radiate  from  the  network  of  branching  channels.  These  delta  lobes 
are  evolved  from  shallow  sandbars  that  rose  above  mean  sea  level  and 
emerged  as  subaerial  land. 

173.  Adams  and  Baumann  (1980)  indicated  that  in  Atchafalaya  Bay 
the  central  area  between  two  deltas  is  apparently  broad  and  deep 
enough  to  transport  the  riverflow  it  receives  without  forming  a 
discrete  channel.  This  fact  suggests  that  in  our  study,  the  analyt¬ 
ical  results  derived  from  numerical  procedures  for  the  Lower 
Atchafalaya  River  Outlet  and  the  Wax  Lake  Outlet  separately  can  be 
combined  linearly  to  represent  the  total  delta  growth  in  the  bay. 

This  is  similar  to  Letter's  (1982)  approach,  in  which  the  bay  area  was 
roughly  divided  into  two  areas,  one  for  each  of  the  two  outlets. 

Figure  33  shows  an  estimation  of  the  total  subaerial  land  in 
Atchafalaya  Bay  derived  from  our  analytical  approach.  Figure  34 
presents  the  prediction  of  volume  extent  of  Atchafalaya  River  delta 
growth. 
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Figure  33.  Analytical  prediction  of  subaerial  land  in  Atchafalaya  Bay  (Case  4) 


PART  VII:  SENSITIVITY  ANALYSIS  AND  RESULTS 


Sensitivity  Analysis  of  Various  Parameters 


174.  The  analytical  solutions  of  jet  characteristics  derived  in 
PART  V  are  based  on  the  theory  of  turbulent  plane  jet  under  a  number 
of  assumptions.  The  numerical  approaches  to  predict  the  delta  growth 
presented  in  PART  VI  depend  upon  a  number  of  estimated  parameters.  It 
is  thus  essential  to  study  the  importance  of  various  physical  param¬ 
eters  in  the  problem  formulation  and  to  show  the  balance  between  the 
numerical  procedures  and  the  physical  environment. 

175.  The  sensitivity  analysis  is  used  to  aid  in  the  under¬ 
standing  of  the  dynamics  of  river-delta  interaction;  to  identify  the 
relative  importance  of  various  variables  in  deltaic  processes;  and  to 
test  the  effects  of  the  distributary  network  on  the  outcome  of  pre¬ 
dictions.  The  sensitivity  tests  to  be  conducted  are: 

a.  River  outlet  conditions  (h^,  b^,  u^) . 

b.  Sediment  concentration  and  settling  velocity  (c^,  . 

£.  Bay  bottom  slope,  friction,  and  lateral  entrainment 
(a,  f,  e). 

176.  The  sensitivity  analysis  is  performed  in  two  phases.  First, 
each  parameter  is  analyzed  independently  by  holding  other  variables 
constant,  under  the  same  conditions  as  in  the  base  results.  Second, 
several  parameters  are  examined  conjunctively  to  show  their  interrelated 
effects . 

River  outlet  conditions  (h^,  b^,  u^) 

177.  In  this  study,  the  water  depth  of  the  receiving  bay  is  taken 
as  the  same  as  the  water  depth  at  the  river  outlet.  To  examine  the 
influence  of  the  outlet  depth  (h^)  on  the  growth  of  delta,  three 
different  water  depths  of  1.5  m,  1.0  m,  and  0.5  m  for  Case  4  (linear 
depth  with  entrainment,  E  ^  0  and  H  ^  1)  are  considered.  That  is,  in  the 
sensitivity  test,  only  the  outlet  water  depth  (hj^)  in  Case  4  is  changed 
from  2.0  m  to  1.5  m,  1.0  m,  and  0.5  m,  while  all  other  variables  remain 
constant.  The  suspended  sediment  distributions  are  given  in  Plate  7  as 


three-dimensional  plots.  The  results  demonstrate  that  at  shallower  water 
depths,  delta  buildup  is  sharper  in  shape  and  is  limited  to  a  much 
smaller  area  close  to  the  river  outlet. 

178.  The  remaining  outlet  conditions  are  studied  by  observing 
simulated  jet  behavior  under  the  influence  of  width  of  the  outlet  (^q) 
and  outlet  velocity  (Uq) .  A  sensitivity  test  is  made  for  the  half-width 
at  values  of  250  m,  500  m,  750  m,  and  1000  m  while  other  variables  in  Case 
4  remain  constant.  The  results  of  the  effects  on  jet  width,  jet 
velocity,  and  sediment  concentration  are  plotted  in  Plate  8.  It  is 
concluded  that  the  larger  the  river  outlet,  the  faster  the  jet  spreads 
laterally  and  the  faster  the  current  velocity  and  sediment  concentration 
diminish  longitudinally.  Also,  from  Plate  9  it  is  shown  that  as  the 
outlet  velocity  (u^)  exceeds  0.5  m/sec,  the  sediment  concentration  (c)  is 
not  substantially  influenced  by  higher  values  of  u^  (u^  =  1.0,  1.5 
m/sec) . 

Sediment  concentration  and  settling  velocity  (c^,  w^) 

179.  A  continous  point  source  of  sediment  issuing  from  the  river 
outlet  is  assumed  in  this  study.  The  supply  of  sediment  to  Atchafalaya 
Bay  has  been  changing  both  in  volume  and  character  over  the  past  decade. 

In  a  regressional  analysis,  Letter  (1982),  based  on  the  50-year  extrap¬ 
olation  hydrograph  at  Siramesport,  computed  the  maximum  and  minimum 
sediment  yield  (113  and  38  million  tons/year)  corresponding  to  the 
maximum  and  minimum  discharges  (310,000  and  139,000  cfs)  and  obtained 
maximum  and  minimum  sediment  concentrations  of  365  ppm  and  275  ppm. 

180.  A  sensitivity  test  is  conducted  for  the  mean  concentration 
(Cq)  (Cq  =  100,  200,  300,  400,  and  600  ppm)  in  Case  4.  The  three- 
dimensional  plots  of  sediment  distributions  are  given  in  Plate  10.  The 
deposited  sediments  become  thicker  as  the  sediment  concentrations  are 
increased;  however,  delta  growth  in  all  instances  is  restricted  to  the 
proximity  area  of  the  outlet. 

181.  The  settling  velocity  of  suspended  sediment  (Wq)  is  in  part  a 
function  of  particle  size  (Figure  23);  the  situation  is  further  compli¬ 
cated  by  the  aggregation  of  the  suspended  cohesive  materials  (Van  Heerden, 
Wells,  and  Roberts  1981).  The  influence  of  the  settling  velocity  is 


103 


investigated  for  =  0.01,  0.05,  0.1,  0.5,  and  1.0  mm/sec.  The  results 
for  the  dimensional  center-line  concentration  are  given  in  Plate  11.  The 
dependence  of  on  c  is  shovm.  The  sediment  deposition  pattern  for  set¬ 
tling  velocities  of  0.01,  0.1,  and  1.0  mm/sec  are  presented  in  Plate  12. 
It  demonstrates  that  smaller  settling  velocities  result  in  more  uniform 
deposition  patterns,  and  that  larger  settling  velocities  produce  patterns 
sharper  in  shape  and  limited  to  a  smaller  area  close  to  the  river 
outlet. 

Bay  bottom  slope,  friction, and  lateral  entrainment  (a,  f,  e) 

182.  In  Atchafalaya  Bay,  Adams  and  Baumann  (1980)  estimated  that 
the  bay  has  a  slight  slope  of  0.00015,  approximately  0.8  ft  per  mile.  In 
this  study,  a  sensitivity  test  is  made  using  slope  values  of  0.00001, 
0.00005,  0.0001,  0.0002,  and  0.001  with  other  parameters  in  Case  4 
remaining  constant.  The  results  of  the  effects  on  jet  width,  jet 
velocity,  and  sediment  concentration  are  shown  in  Plate  13.  The  plots 
demonstrate  that  the  influence  of  the  bottom  slope  is  minor  except  for  a 
value  of  0.001. 


183.  The  dimensional  jet  width  (B),  the  center-line  velocity  (U), 
and  sediment  concentration  (C)  are  displayed  in  Plate  14  for  Case  1  (E=0, 
H=l),  Case  2  (E=0,  Hj^l),  Case  3  {EfO,  H=l),  and  Case  4  (E/0 ,  H/^l)  .  It  is 
seen  that  increasing  the  bottom  slope  causes  a  narrower  and  elongated 
riverine  jet  (Case  1  vs  Case  2;  Case  3  vs  Case  4),  and  that  the  effect 

of  bottom  slope  counteracts  the  effect  of  the  lateral  entrainment  (Case  1 
vs  Case  4) .  The  velocity,  however,  does  not  exhibit  any  significant 
variation  among  the  cases.  The  sediment  concentration  is  affected  both 
by  the  bottom  slope  and  the  lateral  entrainment. 

184.  The  analytical  solutions  of  jet  width,  center-line  velocity, 
and  sediment  concentration  all  show  a  strong  dependence  on  the  bottom 
friction  (f).  In  fact,  the  jet  width  grows  and  jet  velocity  decays 
exponentially  along  the  longitudinal  distance.  The  values  of  the  Darcy- 
Weisbach's  coefficient  are  taken  as  0.001,  0.002,  0.003,  0.004,  and  0.006 
in  the  sensitivity  test  while  all  other  parameters  remain  the  same  as  in 
Case  4.  The  plots  of  jet  width,  center-line  velocity,  and  concentration 
are  shown  in  Plate  15.  It  is  inferred  from  these  figures  that  the  bottom 


friction  plays  an  important  role  in  the  jet  dynamics.  V/hen  the  friction 
is  larger,  the  jet  expands  laterally  much  faster  as  it  faces  the  bottom 
resistance  and  loses  its  momentum  much  more  quickly.  It  also  indicates  a 
rapid  decrease  in  the  sediment  concentration  with  increasing  values  for 
the  bottom  friction. 

185.  Ozsoy  (1977)  stressed  the  importance  of  the  effects  of  lateral 
entrainment  on  tidal  inlet  characteristics.  To  study  the  effects  of 
entrainment  on  the  characteristics  of  the  river  delta,  four  values  of  the 
entrainment  coefficient  e  =  0.0375,  0.075,  0.150,  and  0.300  are  examined. 
Results  are  shown  in  Plate  16  for  jet  width,  center-line  velocity,  and 
concentration,  respectively.  It  is  seen  that  for  a  river-delta  system 
the  entrainment  mechanism  also  plays  a  role. 


The  Relative  Role  of  Physical  Parameters 


186.  The  sensitivity  analyses  conducted  in  previous  sections  shed 
some  light  on  the  relative  importance  of  various  parameters  in  the  study 
of  river-delta  interaction.  Their  orders  of  apparent  importance  are 
listed  below: 

a.  The  bottom  friction  (f)  influences  jet  flow  to  the 
greatest  extent.  The  jet  loses  its  momentum  due  to 
high  friction,  and  expands  its  width  at  a  faster 
rate  than  in  the  case  of  low  friction. 

b.  The  settling  velocity  (w^)  affects  the  transport  of 
suspended  sediment  to  a  large  extent.  With  higher 
settling  velocity,  the  center-line  sediment  concentration 
drops  more  rapidly  and  the  delta  lobe  is  smaller  and 
closer  to  the  river  outlet. 

£.  The  width  of  the  outlet  (b„),  as  considered  in 

simulating  the  process  of  branching  channels,  has 
significant  influence  on  the  shape  of  subaerial  land. 

The  wider  the  outlet,  the  more  rapidly  current  velocity 
diminishes  with  increasing  longitudinal  distance,  which 
results  in  much  faster  lateral  spreading  of  the  jet,  a 
natural  cause  of  delta  buildup. 

d.  The  lateral  entrainment  (e)  affects  the  jet  charac¬ 
teristics.  Jet  width  increases  with  increasing 
entrainment.  Consequently,  the  jet  velocity  and 
sediment  concentration  decrease  as  the  jet  width 
expands . 
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e.  The  influence  of  the  outlet  depth  (h-)  on  jet 
characteristics  can  be  examined  by  utilizing  the 
aspect  ratio  (Jirka  1981),  defined  as  the  ratio 
depth  to  width,  "hQ/bQ**.  Higher  values  of  this  ratio, 
increasing  h^  or  decreasing  b^,  result  in  a  much 
narrower  and  more  elongated  delta  form. 

f .  Increasing  the  center-line  velocity  (u  )  and  the  sediment 
concentration  (c^)  will  develop  delta  lobes  with  thicker 
and  more  elongated  shapes. 


Growth  Curves  of  Subaerial  Land  in  Atchafalaya  Bay 

187.  In  this  study,  the  values  used  in  the  base  runs  are  mean 
values  of  various  physical  parameters  taken  from  the  available  literature 
(PARTS  II  and  VI).  These  quantities  may  vary  throughout  the  period  of 
predicted  growth  depending  on  the  nature  of  the  system  under  considera¬ 
tion.  To  circumvent  this  situation,  the  most  sensitive  parameters  are 
used  to  simulate  the  greatest  and  least  potential  growth  of  subaerial 
land. 

Simulated  growth  of  Atchafalaya  River  deltas 

188.  The  physical  quantities  that  are  of  major  importance  to  the 
Atchafalaya  River-Bay  system  are  the  coefficient  of  bottom  friction  (f) , 
the  settling  velocity  of  sediment  particles  (Wq),  and  the  sediment 
concentration  (Cq) .  Since  high  friction  causes  the  jet  to  spread 
laterally,  high  settling  velocity  results  in  smaller  delta  lobes,  and  low 
sediment  concentration  develops  a  thinner  subdelta,  the  environment  for 
least  growth  can  be  simulated  by  increasing  f  and  w^  and  decreasing  the 
Cq  value.  In  contrast,  the  high  growth  environment  can  be  best  simulated 
by  low  f,  low  Wq,  and  high  c^  values. 

189.  The  following  environments  are  simulated  for  the  growth  of 
subaerial  land  in  Atchafalaya  Bay: 

a.  Base  results 

f  =  0.001 

Wq  =  0.05  ram/sec 

Cq  =  300  ppm 
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b.  Slow  growth 


f  =  0.004 


Wq  =  0.1  mm/sec 


Cq  =  200  ppm 


c.  Fast  growth 


f  =  0.001 


Wq  =  0.03  mm/sec 


Cq  =  400  ppm 


The  results  of  the  simulated  growth  patterns  are  summarized  in  Tables  7, 
8,  9,  and  10  and  plotted  in  Plate  17. 


Comparison  with  other  predictions 


190.  The  growth  of  the  Atchafalaya  River  deltas  has  been  predicted 
by  various  investigators.  Shlemon  (1972),  based  on  sediment  measurements 


made  in  the  outlets,  predicted  a  growth  rate  of  5.5  to  7.5  mi^/year  (14 


2  2 
to  18  km  /year),  with  the  delta  covering  an  area  of  290  to  350  mi  (750 


to  900  km  )  by  the  year  2020.  His  prediction  was  referred  to  a  -3  ft 
mean  low  gulf  contour. 

191.  Adams  and  Baumann  (1980),  following  the  same  empirical 
approach  as  Shlemon,  estimated  that  the  Atchafalaya  Bay  will  be  filled  to 
an  average  depth  of  2  ft  below  mean  sea  level  in  a  time  period  of  approx¬ 
imately  40  years. 

192.  Wells,  Chinburg,  and  Coleman  (1984),  based  on  the  generic 


analysis  of  existing  deltas,  projected  that  by  the  year  2030  a  new  sub- 

2 

aerial  land  mass  will  be  created  in  the  bay  ranging  from  59  to  132  mi 


2  2  2 

(150  to  337  km  )  with  81  mi  (208  km  )  representing  the  expected  land  in 


50  years  under  average  flood  conditions.  They  estimated  a  growth  rate  of 


2  2 

1.6  mi  /year  (4.0  km  /year);  their  study  was  referred  to  the  mean  sea  level. 


193.  Letter  (1982)  developed  a  regressional  model  and  predicted 
that  within  50  years  (by  the  year  2027)  the  delta  will  evolve  gulfward  of 


Eugene  Island.  The  total  volume  of  the  deposited  sediment  is  estimated 


at  58  billion  ft  ;  and  the  delta  mass  volume  (based  on  the  -3  ft  NGVD)  is 
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about  17.6  billion  ft' 


194.  The  total  subaerial  land  development  predicted  analytically 
in  this  study  is  displayed  in  Figure  35  together  with  the  prediction 
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Figure  35.  Prediction  of  subaerial  land  in  Atchafalaya  Bay 

curves  developed  by  Shlemon  (1972)  and  Wells,  Chinburg,  and  Coleman  (1984) 
A  contour-type  map  for  approximate  delta  front  advancement  is  depicted  in 
Figure  36.  Our  prediction,  based  on  analytical  results,  shows  an  average 
growth  rate  of  7.6  km  /yr,  ranging  from  6.0  to  10.0  km  /yr  as  simulated 
for  the  slow  growth  and  fast  growth  environments,  respectively.  On  a 
volume  basis,  the  average  volume  of  sediment  deposition  is  predicted  at 
16  X  10^  m^/yr,  with  a  range  of  12  x  10^  m^/yr  to  23  x  10^  m^/yr. 


PART  VIII:  SUMMARY,  LIMITATIONS,  AND  CONCLUSIONS 


Summary 

195.  An  analysis  has  been  made  to  aid  in  understanding  the  various 
phenomena  associated  with  turbulent  plane  jets  issuing  from  river  outlets 
and  discharging  into  a  quiescent  receiving  bay.  An  integrated  form  of 
the  hydrodynamic  equations  of  flow  continuity  and  momentum  balance, 
coupled  with  the  mass  transport  equation,  has  been  formulated  into  a 
two-dimensional  spatial  and  quasi-steady  temporal  domain. 

196.  A  similarity  function  in  the  form  of  exponential  and  poly¬ 
nomial  expressions  was  chosen  for  the  velocity  and  sediment  concentration 
profiles.  The  lateral  entrainment  is  expressed  as  a  function  of  the  jet 
center-line  velocity,  A  closed-form  analytical  solution  is  obtained  in 
terms  of  the  normalized  dependent  variables  of  the  jet  width,  the  center- 
line  jet  velocity,  and  the  center-line  sediment  concentration.  The 
solutions  for  cases  of  constant  depth  and  of  linearly  varying  depth,  both 
with  and  without  entrainment,  are  presented. 

197.  From  these  normalized  solutions,  the  thickness  of  deposited 
sediment  is  calculated  under  quasi-steady  state  conditions.  The  Statis¬ 
tical  Analysis  System  of  computer  graphics  (SAS/GRAPH  1981)  is  used  to 
display  the  three-dimensional  sediment  deposition  patterns.  These 
sediment  patterns,  together  with  a  stepwise  numerical  procedure  simu¬ 
lating  the  process  of  branching  channels,  are  then  used  to  estimate  the 
areal  extent  and  the  volume  of  deposition  for  the  Atchafalaya  deltas 
under  a  variety  of  conditions.  Sensitivity  analyses  are  performed  to 
assess  the  relative  importance  of  various  parameters  in  the  river-delta 
system. 


Limitations 

198.  There  are  certain  limitations  inherent  in  the  application  of 
the  analytical  model: 


a.  The  governing  equations  for  flow  and  mass  transport  are 
developed  for  a  shallow-water  environment  in  the  absence 
of  tidal  and  wave  currents.  Therefore  the  simplified 
equations  for  turbulent  jets  will  not  be  representative 
of  a  regime  subjected  to  tidal  and  wave  action. 

b.  The  similarity  profiles  for  the  velocity  and  sediment 
concentration  must  be  used  in  order  to  obtain  a  set  of 
closed-form  analytical  solutions. 

£.  The  deposition  function  used  for  sediment  dispersion  and 
settling  in  shallow  water  does  not  take  the  erosion 
process  into  consideration.  The  analytical  approach 
is  not  capable  of  addressing  the  problem  of  the 
resuspension  of  sediment  in  the  bay  and  the  reworking 
of  delta  deposits  by  physical  processes  offshore. 

d.  In  studying  the  river  delta  development,  the  sediment 
deposition  is  assumed  to  be  quasi-steady  state;  the 
thickness  of  the  deposited  sediment  layer  is  assumed  to 
vary  linearly  with  time. 

e.  The  bifurcation  processes  are  simulated  in  the  form  of 
turbulent  jets  at  the  new  channel  outlets;  at  each  level 
of  bifurcation,  an  equal  amount  of  flow  and  sediments  is 
assumed  to  be  distributed  to  the  branching  channels. 

f .  The  values  of  various  parameters  are  based  on  time- 
averaged  quantities.  The  analysis  is  projected  to  yield 
a  gross  estimate  of  areal  and  volume  extent  of  the 
Atchafalaya  River  delta. 

g.  The  analytical  results  are  generated  from  local  data 
within  the  bay.  Predicting  the  delta  development  beyoi.d 
the  bay  will  be  less  accurate,  since  the  governing 
equations  of  turbulent  plane  jets  derived  for  shallow 
receiving  waters  will  not  be  applicable  to  the  deeper 
water  offshore. 


Conclusions 


199.  The  primary  sources  of  energy  for  the  development  of  the 
Atchafalaya  River's  deltas  in  the  bay  are  the  natural  resources  of  the 
river,  that  is,  the  suspended  sediments  that  form  the  delta  and  the  river 
discharge  that  carries  the  sediment  load.  The  river-delta  interaction  is 
a  complicated  phenomenon.  The  analytical  work  reported  herein  is  a 
simplified  representation  of  complexity  in  predicting  the  evolution  of 
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the  Atchafalaya  River  delta.  It  covers  the  essential  features  of  natural 
processes  leading  to  the  growth  of  the  delta. 

200.  The  following  conclusions  are  drawn  from  this  study: 

a.  The  depth  of  the  river  outlet  and  the  bottom  friction  in 
the  bay  play  an  important  role  in  determining  the 
spreading  of  river  effluents. 

b.  The  process  of  channel  bifurcation  has  significant 
influence  in  the  shape  and  the  area  of  subaerial  land. 

c.  The  settling  velocity  of  sediment  particles  has  the  most 
impact  on  the  volume  of  delta  lobes. 

d.  The  center-line  velocity  and  sediment  concentration 
control  both  the  area  extent  and  volume  deposition 
of  river  deltas. 

e.  The  lateral  entrainment  also  plays  a  role  in  the 
development  of  river  deltas  in  the  shallow  and  wide 
Atchafalaya  Bay. 

f.  The  predictions  of  the  future  growth  of  the  Atchafalaya 
River  deltas  in  this  study  are: 

(1)  The  Lower  Atchafalaya  River  Outlet  Delta  will 
grow  approximately  5.0  km^/yr;  the  Wax  Lake  Outlet 
Delta  will  grow  at  a  rate  of  2.6  km^/year. 

(2)  The  total  growth  of  subaerial  land  of  the 
Atchafalaya  River  deltas  is  expected  to  be 
7.6  km^/year. 

g.  It  will  take  about  50  years  for  the  Atchafalaya  River 
deltas  to  reach  the  Point  Au  Fer  Shell  Reef. 

h.  The  results  of  delta  growth  predicted  in  this  analyt¬ 
ical  study  are  commensurate  with  the  predictions  made 
by  others. 
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A  three-branch  distribution  pattern  Is  used  for  the  second  and  third  levels  of  bifurcation; 
for  the  remaining  levels,  the  distributary  networks  are  characterized  by  two-branch  channel 
patterns. 
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APPENDIX  A: 


GAUSSIAN  QUADRATURE  NUMERICAL  INTEGRATION 
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1.  To  evaluate  the  integral  of  a  function,  f(x),  over  a  finite 
interval  (-1,1),  Gaussian  quadrature  numerical  integration  (Beyer  1978)* 
is  used.  The  Gaussian  quadrature  formula  has  the  form 

1 


I 


n 

f(x)dx  =  1  (H.f(x.)] 
i=i 


(Al) 


where  are  the  weights;  the  abscissas,  x^,  occur  in  pairs  symmetrically 

placed  with  respect  to  the  origin. 

2.  For  an  eight-point  numerical  integration  (n  =  8),  the  values 

for  X.  and  H.  are  given  as  follows: 

1  1 


±  0.9602899 
±  0.7966665 
±  0.5255324 
±  0.1834346 


H. 

1 

• 

0.1012285 

0.2223810 

0.3137066 

0.3626838 


3.  A  simple  FORTRAN  coding  for  the  Gaussian  quadrature  formula  is 
listed  in  Table  Al.  The  values  of  the  following  integrals,  cited  in  PART 
V,  are  found: 


J  (1  -  s^)  exp  (-s^)  ds  =  1.1147022 


-1 


/  (1  -  s^)  ^  exp  (-  |s^)  ds  =  0.9494728 


-1 


J  (1  -  s^)^  exp  (-  |s^)  ds  =  1.3974447 


-1 


/  (1  -  s^)^  exp  (-  |s^)  ds  =  0.7597358 


-1 


*References  cited  in  the  appendices  are  included  with  those  for  the  main 
body  of  the  report,  starting  on  page  113. 


Al 


Eight- Point  Numerical  Integration 


GAUSSIAN  QUADRATURE  EIGHT  POINT  NUMERICAL  INTEGRATION 
REFERENCE:  BEYER,  W.H..  1978.  "HANDBOOK  OF  MATHEMATICAL 

SCIENCE".  5TH  EDITION.  CRC  PRESS. 

GAUSSIAN  QUADRATURE  FORMULA: 

WHERE  X  =  ABSCISSAS  AND  H  =  WEIGHTS 

DIMENSION  X(8) ,H(8) .X1(8) ,H1 (8) 

DATA  FI  ,F2,F3,F4/4i«0.0/ 

REAL*8  X 1 /O . 9602899 . 0 . 7966665 . 0 . 5255324 ,0.1 834346 , 

-0 . 1 834346 . -0 . 5255342 . -0 . 7966665 , -0 . 9602899/ . 

H 1 /O . 1 0 1 2285 . 0 . 22238 10,0.31 37066 , 0 . 3626838 . 

0 . 3626838 ,0.31 37066 . 0 . 22238 10,0.101 2285/ 

DO  9  1=1 ,8 
X(I)=X1 (I) 

9  H(I)=H1(I) 

WRITE(6. 1 ) 

1  FORMAT (5X. 'THE  FOLLOWING  INTEGRALS  ARE  COMPUTED:') 

DO  11  1=1,8 

11  F1=F1+H(I) *( 1 .-X(I)»X(I) )*EXP(-X(I)wX(I)  ) 

WRITE(6, 1 0)F1 

DO  12  1=1,8 

12  F2=F24-H(I)»<(  1  .-X(I)*X(I)  )**  1 .5*EXP(-1  .5*X(I )  *X(I)  ) 
WRITE(6,20)F2 

DO  13  1=1,8 

13  F3=F3+H(I)'«(  1  .-X(I)*X(I)  )*»0.5*EXP(-0.5*«X(I)*X(I)  ) 
WRITE(6,30)F3 

DO  14  1=1,8 

14  F4=F4+H(I)*( 1 .-X(I)«X(I) ) * *2 . 5*EXP ( -2 . 5*X ( I) *X ( I )  ) 
WRITE(6.40)F4 

STOP 

10  FORMAT (  10X, 'K  1  )  =',F10.7) 

20  FORMATC 10X, ’1(3/2)  =',F10.7) 

30  FORMATC 1  OX,  ’ I ( 1/2)  =',F10.7) 

40  FORMATC 1  OX,  ’1(5/2)  =',F10.7) 

END 

THE  FOLLOWING  INTEGRALS  ARE  COMPUTED: 

1(1)  =  1.1 147020 

1(3/2)  =  0.9494730 
1(1/2)  =  1.3974440 
1(5/2)  =  0.7597359 


APPENDIX  B:  SOLUTION  FOR  BERNOULLI -TYPE  EQUATIONS 


1.  This  appendix  shows  the  procedure  for  reducing  a  Bernoulli- type 
equation  into  a  linear  first-order  differential  equation.  Bernoulli-type 
equations  have  the  form 


^  +  p(x)y  =  q(x)y" 


This  basic  equation  is  modified  by  setting 


2  =  y 


(B2a) 


which  transforms  into  the  equation 


y  =  2 


(B2b) 


The  first  term  of  Equation  Bl  can  then  be  expressed  as 


^  ^  ^  _  1  1-n  ^ 

dx  d2  dx  1-n  dx 


Substituting  Equations  B2a  and  B3  into  Equation  Bl  and  simplifying, 
we  have 


^  +  (1  -  n)p(x)2 


=  (1  -  n)q(x) 


(B^a) 


^  +  P(x)2  =  Q(x) 


(B4b) 


l-v- 


vv-; 


m 


C<' 


2.  Equation  B4b  is  the  general  form  of  a  linear  first-order 
differential  equation.  By  the  method  of  integrating  factors  (Wylie 
1951),  the  general  solution  of  Equation  B4b  is 


|jQ(3c)dx  +  Cj] 


(B5) 


where 


p  =  integrating  factor 
=  exp  f"p(x)dx 


(B6) 


and 


C  =  constant  of  integration 


3.  Finally,  using  the  relationship  of  Equation  B2b,  the  standard 
solution  for  the  Bernoulli-type  equation  is  obtained: 


y  =  2 


1 

1-n 


■/' 


.-1 


y  =  { [exp  j  P(x)dx]  [  exp  /  P(x)dx  Q(x)dx  +  C] } 


.  1 

1-n 


(B7) 


where 


P(x)  =  (1  -  n)p(x) 
Q(x)  =  (1  -  n)q(x) 


(B8a) 

(B8b) 


The  Solution  of  Equation  53  in  PART  V 


In  PART  V,  Equation  53  is  a  Bernoulli-type  ordinary  differential 
equation : 


dU 


dr 


+  FU  =  -  A  exp  (Fr)  U' 


(B9) 


B2 


The  solution  of  the  above  equation  can  be  obtained  by  equating 


y  =  U 


n  =  3 


P(x)  =  (1  -  n)p(x)  =  -2F 

Q(x)  =  (1  -  n)q(x)  =  2A  exp  (Fr) 


and  computing  the  integrating  factor 


p  =  expj  P(x)dx 

=  exp  J"  -2Fdr 

=  exp  (-2Fr) 


U  =  lexp  (2Fr)  i-  ^  exp  (-Fr)  C]}  ^ 


Utilizing  the  initial  condition 


"lr=0=  ‘  = 


we  obtain 


C  =  1  + 


Finally,  the  solution  of  the  Bernoulli  equation  is 


U  =  {exp  (2Fr)(-^  exp  (-Fr)  +  1  +  ^1}’^ 


U  =  [-  ^  exp(Fr)  +  (1  +  -j:)  exp  (2Fr)l 


(B16J 


This  is  the  jet  velocity  ia  nondimensional  form  for  the  case  of  constant 
depth.  It  is  cited  as  Equation  54  in  PART  V. 

The  Solution  of  Equation  72  in  PART  V 

Equation  72,  appearing  in  PART  V,  is  also  a  Bernoulli-type 
ordinary  differential  equation; 


+  ar^(D  -  1)H‘^  (HU)  = 
dr  hg 


(B17) 


The  solution  is  obtained  as  in  the  previous  example,  using 


y  3  HU 


X  3  L- 


n  =  3 


P(x)  =  (1  -  n)p(x)  = 


‘"O  -1 

-2ar^(D  -  1)H  ' 


Q(x)  =  (1  -  n)q(x)  = 


The  integrating  factor  is  computed  as 


(Bl8a) 


(Bl8b) 


tj  =  exp  /  P(x)dx 


=  exp 


y*  P(x)d: 


(D  -  l)H"^dr 


(B19a) 


(B19b) 


Note  that 


H  = 


1  +  ar^r 


and 


dH  = 


V 

V' 


Thus 


M  =  exp 


f-ia-  i)f 


and 


HU  = 


Utilizing  the  initial  conditions 

HU 


=  1  =  ^  ♦  c 

r=0  2  -  D  abg  ^ 


the  constant  of  integration  is  found  to  be 


C  =  1  - 


2A  ^0 


2  -  D  ab. 


The  final  solution  of  this  Bernoulli  equation  is 


HU 


_  r„2(D-l).  2A  „2-D  .  ,  2A  ‘'O 

-  iH  ^  “  +  1  -  2-t-d  ^ 


or 


HU  = 


(2A)-Si'-'>ij§^  ■  '>  "  k''‘‘ 


The  above  equation  appears  as  Equation  73  in  PART  V. 


APPENDIX  C:  COMPUTER  GRAPHICS 


I  1.  Behavior  of  variables  used  in  the  system  of  equations  employed 

I 

>  to  describe  sedimentation  patterns  in  Atchafalaya  Bay  is  Illustrated  in 

the  report  by  the  use  of  Statistical  Analysis  System  plotting  subroutines 
(SAS/GRAPH  1981).  Two-dimensional  plots  are  used  to  display  changes  in 
jet  center-line  values  of  selected  variables  as  distance  from  the  outlet 
increases;  three-dimensional  plots  are  used  to  depict  values  of  the 
variables  over  the  entire  jet. 

2.  The  GPLOT  procedure,  a  subprogram  of  SAS/GRAPH,  was  used  to 
produce  the  majority  of  the  two-dimensional  plots.  Variables,  entered  as 
X  and  y  coordinate  pairs,  are  plotted  exactly;  the  plotting  routine  in¬ 
cludes  an  interpolation  feature  (SPLINE)  that  produces  a  smooth  curve  fit 

I 

to  the  data  values. 

I  3.  A  second  SAS/GRAPH  subprogram,  the  G3D  procedure,  plots  the 

I  values  of  three  variables  and  generates  a  three-dimensional  surface.  The 

I  three  variables  used  to  illustrate  deposition  patterns  are  jet  velocity, 

I  sediment  concentration,  and  sediment  thickness.  Each  of  these  variables 

I 

'  is  plotted  versus  normalized  longitudinal  (r)  and  lateral  (s)  coordinates 

by  the  program,  which  again  interpolates  to  produce  a  smooth  surface  fit 
to  the  input  data. 


Computer  Program  Cl 

4.  This  computer  program  contains  plotting  instructions  for  the 
generation  of  a  two-dimensional  center-line  sediment  thickness  profile 
using  the  GPLOT  procedure.  The  data  used  to  create  the  plot  were  taken 
from  the  Case  4  conditions  for  delta  development  (E  7^  0  and  linear  H) , 
after  two  years  of  deposit,  with  varying  sediment  concentrations.  Pre¬ 
dicted  sediment  thickness  at  forty  equidistant  locations  along  the  center 
line  of  the  jet  was  entered  with  corresponding  longitudinal  distance 
from  the  outlet;  data  points  are  indicated  on  the  plot  by  designated 
symbols  representing  different  sediment  concentrations.  The  resulting 
plot  is  shown  in  Figure  Cl. 


Cl 


COMPUTER  PROGRAM  Cl 


/•/STEPl  EXEC  SAS 
//SYSIN  DD  « 

DATA  ONE; 

U0=1;  W0=0.05E-3;  F=0.001;  80=500; 

H0=2.0;  A=0.000l;  AIt=0.948;  AI2=1.397;  AI3=0.76; 

TIME=2. 0*365. *86400. ;  DS=0.4; 

ARRAY  ROOT  ( I )  ROOT  1  -  ROOT  1 1 ; 

ARRAY  WT  (I)  WTl  -  WTII* 

ROOT 1 =0 . 00533697 ;  Ro6t2=0 . 04735263 ;  R00T3=0 . 1 2782 1 06 ; 
ROOT4=0 . 2399 1861;  ROOT5=0 . 374 1 3950 ;  ROOT6=0 .51910186; 
ROOT7=0. 6625 1286;  R00T8=0. 79221 088;  ROOT9=0 . 897 1 9596 ; 

ROOT  1 0=0 . 96855602 ;  ROOT 11  =  1. 00000000 ; 

UT 1=0. 29 169804;  UT2=0 . 28548098;  UT3=0 . 273 1 7938 ; 

UT4=0 . 25505539 ;  UT5=0 . 23 1 49529 ;  WT6=0 . 20300 115; 

WT7=0. 17018012;  WT8=0 . 1 3373 1 2 1 ;  MT9=0 . 09442893; 

WT 1 0=0 . 05309 150;  UT 11 =0 . 0086580 1 ; 

UCR  =  UO; 

D  =  F/A; 

U  =  BO*WO/UO/HO; 

D  =  F/A;  E=0.075; 

AA  =  1 .794*E; 

ARRAV  CO(I)  C01  C05; 

C0 1  =  1 00 . ;  C02=200 . ;  C03=300 . ;  C04=400 . ;  C05=600 . ; 

ARRAY  CA(I)  CA1  -  CA5; 

ARRAY  CD(I)  CD1  -  CD5; 

ARRAY  UD(I)  UD1  -  UD5; 

ARRAY  H  (I)  HI  -  HI  1 ; 

ARRAY  CFU  (I)  CFU1  -  CFU1 1 ; 

ARRAY  U  (I)  U1  -  U1 1 ; 

ARRAY  SA  (I)  SA1  -  SA1 1 ; 

ARRAY  SB  (I)  SB1  -  SB1 1 ; 

DO  R  =  0.0  TO  20.0  BY  0.50; 

HH  =  1.  +  A*BO/HO*R; 

CF  =  H0/A/B0/(2.-D)*(HH**<2.-D)  -  1.)  +  0.5/AA; 

UU  =  1 ./SQRT(2. *AA)  /HH**D  /  SQRT(CF) ; 

BB  =  2. *AA*HH**(D-1 )  *  CF; 

DO  OVER  ROOT;  H  =  1.  +  A*BO/HO*R*SQRT (ROOT) ;  END; 

DO  OVER  H;  CFU=  H0/A/B0/(2. -D) * (H»* (2 . -D)  -  1.)  +  0.5/AA;  END 
DO  I  =  1  TO  11; 

U  =1./SQRT(  2.«AA)  /  H**D/SQRT(CFU) ; 

END; 

G1=  -AI2*W/AI 1 »R/2. ; 

G2=  AI3-W/AI 1/1 .0/1 .0*R/2. ; 

DO  I  =  1  TO  11; 

SA  =  WT/H/U  ; 

SB  =  WT*U/H  ; 

END; 

SUM1  =  SUMCOF  SA1  -  SAID; 

SUM2  =  SUMCOF  SB  1  -  SB  1  1)  ;' 


C2 


COMPUTER  PROGRAM  Cl  (CONTINUED) 


C  =  EXP(  G1*SUM1  +  G2»SUM2  -  LOG(HH*UU*BB)  ); 

DO  I  =  1  TO  5; 

CA  =  CO  «  C; 

UD=UO»UU; 

CD  =  WO-IOO.^CA*!  .0E-6»(  1  .  -  UD>«UD/UCR/UCR )  *TIME ; 

END; 

TA  =  CD1/DS; 

T3  =  CD2/DS; 

TC  =  CD3/DS; 

TD  =  CD4/DS; 

TE  =  CD5/DS; 

OUTPUT; END; 

PROC  GPLOT;  PLOT  TAxR  TB*R  TC«R  TD-R  TE*R/OVERLAY; 

SYMBOL!  I=SPLINE  V=DIAMOND; 

SYMB0L2  I=SPLINE  V=TRI ANGLE; 

SYMBOLS  I=SPLINE  V=STAR; 

SYMB0L4  I=SPLINE  V=PLUS; 

SYMBOLS  I=SPLINE  V=HASH; 

TITLE  1  ; 

TITLE2  ; 

TITLE3  .F=TRIPLEX  .H=2  CENTERLINE  SEDIMENT  THICKNESS  PROFILE; 
TITLE4  .F=TRIPLEX  .H=l  IN  CM,  TWO  YEARS  DEPOSIT; 

TITLES  .F=TRIPLEX  .H=l  CASE  4:  E  i?  0  AND  LINEAR  H; 

TITLE6  .F=TRIPLEX  .H=1  C0=100:  DIAMOND,  C0=200:  TRIANGLE; 

TITLE7  .F=TRIPLEX  .H=l  C0=300 :  STAR,  C0=400:  PLUS,  C0=600 :  HASH; 
/■ 


CENTER-LINE  SEDIMENT  THICKNESS  PROFILE 

IN  CM,  TWO  YEARS  DEPOSIT 
CASE  4:  E  0  AND  LINEAR  H 
C0=100:  DIAMOND.  C0=200:  TRIANGLE 
C0=300:  STAR,  C0=400:  PLUS,  C0=600:  HASH 


5.  Computer  program  C2  contains  plotting  instructions  for  genera¬ 
tion  of  a  three-dimensional  representation  of  sediment  thickness  using 
the  G3D  procedure;  the  resulting  plot  is  shown  in  Figure  C2.  The  sedi¬ 
ment  deposition  pattern  shown  was  predicted  by  the  closed-form  analytical 
solutions  of  the  two-dimensional  system  of  equations  describing  sediment¬ 
ation  for  Case  4  (E  0  and  linear  H)  after  two  years  of  deposit.  These 
equations  are  used  by  the  G3D  subroutine  to  produce  1651  data  points  for 
generation  of  the  final  three-dimensional  plot. 


COMPUTER  PROGRAM  C2 


//STEPI  EXEC  SAS 
//SYS IN  DD  » 

DATA  ONE; 

U0=1;  W0=0.05E-3;  F=O.OOl;  80=500;  C0=300.; 

H0=2.0;  ^=0.0001 ;  AI1=0.948;  AI2=1 .397;  AI3=0.76; 

TIME=2. 0*365. *86400. ;  DS=0.4; 

ARRAY  ROOT  ( I )  ROOT  1  -  ROOT 1 1 ; 

ARRAY  WT  (I)  WT1  -  WTt 1 ; 

ROOT 1 =0 . 00533697 ;  ROOT2=0 . 04735263 ;  ROOT3=0 . 1 2782 1 06 ; 
ROOT4=0. 23991 86 1 ;  ROOT5=0 . 374 1 3950 ;  ROOT6=0 . 5 1 9 1 0 1 86 ; 
ROOT7=0. 6625 1286;  ROOT8=0 . 7922 1 088 ;  ROOT9=0 . 897 1 9596 ; 

ROOT  1 0=0 . 96855602 ;  ROOT  11  =  1. 00000000 ; 

WT 1=0. 29 169804;  WT2=0 . 28548098 ;  WT3=0 . 273 1 7938 ; 

WT4=0 . 25505539 ;  WT5=0 . 23 1 49529 ;  WT6=0 . 20300 115; 

WT7=0 .17018012;  WT8=0 . 1 3373 121;  WT9=0 . 09442893 ; 

WT 1 0=0 . 05309 150;  WT 1 1 =0 . 0086580 1 ; 

UCR  =  UO; 

D  =  F/A; 

W  =  B0*W0/U0/H0; 

D  =  F/A;  E=0.075; 

AA  =  1 .794*E; 

ARRAY  H  (I)  HI  -  HI  1 ; 

ARRAY  CFU  (I)  CFU 1  -  CFU11; 

ARRAY  U  (I)  U1  -  U1 1 ; 

ARRAY  SA  (I)  SA1  -  SA1 1 ; 

ARRAY  SB  (I)  SB1  -  SB1 1 ; 

DO  R  =  0.0  TO  20.0  BY  0.50; 

DO  S  =  -1  TO  1  BY  0.05; 

HH  =  1.  +  A*BO/HO*R; 

CF  =  H0/A/B0/(2.-D)»(HH*«(2.-D)  -  1.)  +  0.5/AA; 

UU  =  1 ./S0RT(2. *AA)  /HH*»D  /  SQRT(CF) ; 

BB  =  2. •AA*HH»*(D-1 )  ■  CF; 

DO  OVER  ROOT;  H  =  1.  +  A»BO/HO*R»SQRT(ROOT) ;  END; 

DO  OVER  H;  CFU=  H0/A/B0/(2. -D) * ( H»* (2 . -D)  -  1.)  +  0.5/AA;  END; 
DO  I  =  1  TO  11; 

U  =1./SQRT(  2.*AA)  /  H**D/SORT(CFU) ; 

END; 

G1=  -AI2*W/AI 1 »R/2. ; 

G2=  AI3*W/AI1/1 .0/1 .0»R/2. ; 

DO  I  =  1  TO  1 1 ; 

SA  =  WT/H/U  ; 

SB  =  WT*U/H  ; 

END; 

SUMl  =  SUMCOF  SA1  -  SAID; 

SUM2  =  SUMCOF  SBI  -  SBD); 

C  =  EXPC  G1«SUM1  +  G2*SUM2  -  LOG(HH*UU*BB)  ); 

Cl  =  C0*SQRT( 1 .-S*S)/EXP(S*S/2. ) *C; 

CC1=(1.  -  S*S)/EXP(S*S) ; 

UD1=U0»CC1 »UU; 

CD  =  WO* 100. «C1 *  1 .0E-6*( 1 .  -  UD1«UD1/UCR/UCR)*TIME; 

THICK  =  CD/DS; 

OUTPUT; END; END; 


C.6 


COMPUTER  PROGRAM  C2  (CONTINUED) 

TITLE  1  ; 

TITLE2  i 
TITLE3  ; 

TITLE4  ; 

TITLES  ; 

TITLE6  .F=TRIPLEX  . H=2  PLOT  OF  SEDIMENT  THICKNESS: 

TITLE?  .F=TRIPLEX  . H= I  :  IN  CM; 

TITLES  .F=TRIPLEX  . H= 1  CASE  4;  E  =  6  AND  LINEAR  H; 
TITLE9  .F=TRIPLEX  . H= 1  BASE  RESULTS.  2.0  YEARS  DEPOSIT; 
PROC  G3D; 

PLOT  R*<S=THICK; 

// 
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APPENDIX  D;  QUADRATURE  INTEGRATION 


1.  Orthogonal  collocation  is  an  efficient  method  for  solving 
differential  equations  (Villadsen  and  Michelson  1978).  Based  on  this 
method,  the  profile  of  a  dependent  variable  can  be  expressed  in  terms  of 
the  families  of  orthogonal  polynomials  and  their  derivatives  at  colloca¬ 
tion  points.  In  addition  to  its  efficiency  for  solving  differential 
equations  the  integration  of  the  profile  of  dependent  variables  can  be 
closely  approximated  by  a  quadrature  formula  which  gives  a  semianalytical 
solution  of  integration. 


Orthogonal  Polynomial 


2.  The  orthogonal  polynomial  may  be  of  Jacobi,  Legendre,  or 
Chebycheff  type  of  which  the  Jacobi  polynomial  is  found  to  be  the  most 
efficient  (Villadsen  and  Michelson  1978).  The  important  features  of 
orthogonal  polynomials  are  illustrated  in  this  appendix.  The  Jacobi 
family  possesses  the  following  orthogonal  relationship: 


/ 


W(u) 


p  (a,P) 

i 


(u)  du  = 


i  *  j 
i  =  j 


(Dl) 


where  W(u)  is  the  weight  function  for  a  Jacobi  polynomial  with  the  fol¬ 
lowing  expression; 


W(u)  =  (1  -  u)°*  u^ 


(D2) 


in  which  a  and  p  are  constants  greater  that  -1.  The  polynomial  P . (u)  is 

fo  61  ^ 

automatically  fixed  when  the  weight  function  is  given.  P^  (u)  and 

f  o  R 1  ^ 

Pj  (u)  are  the  Jacobi  polynomial  with  degree  i  and  j  respectively. 

3.  A  Jacobi  polynomial  of  degree  n  can  be  expressed  by  the 

following  equation: 


(u)  =  I  (-1)“"^  a  u  (D3) 

i=0  ^ 


where  the  coefficient  a.  can  be  expressed  by  the  following  recurrence 
formula : 


.  _n-i  +  l  n+i  +  a  +  B 
1  i  i  +  P 


(D4) 


with  a^  =  1,  and  i  =  1,2 . n. 

For  n  =  3  and  a=p=0,  Equation  D3  immediately  gives 

p(0,0)  =  20  u^  -  30  u^  +  12  u  -  1  (D5) 

4.  For  convience  of  computer  computation,  Equations  D3  and  D4 
may  also  be  expressed  by  the  following  equations: 


and 


a+p+2  “  ^ 


2  „2 

«  -P 


] 


(2n+a+p-l)‘‘  -  1 

for  n>l 


hi  =  0 


h2  = 


(g+l)  (P+1) 
(a+p+2)^  (a+p+3) 


u  _  (n-1)  (n+g-l)  (n+p-1)  (n+gj-p-l) 

^  >  tor  ii>2 

(2n+g+p-l)  (2n+g+p-2)  (2n+g+p-3) 


(D6) 

(D7) 

(D8) 

(D9) 

(DIO) 


D2 


The  recurrent  evaluation  of  (u)  starts  with  n=l,  (u) 

arbitrarily,  and  (u)=l. 


Quadrature  Integration  Formulae 


5.  In  orthogonal  collocation  the  integration  can  be  conveniently  and 
accurately  performed  by  a  numerical  quadrature  formula.  Two  types  of 
quadrature  formulas  are  described  briefly  in  the  following. 

Gauss-Jacobi  quadrature 

6.  The  Gauss-Jacobi  quadrature  can  integrate  any  polynomial  Y(x) 
up  to  2n-l  degree  by  the  following  equation: 


1  0,  ft  n 

(1  -  x)  xP  Y(x)  dx  =  I  W.  Y(x.) 
'0  i=l  ^  ^ 


/ 


(Dll) 


The  integration  can  be  obtained  closely  by  summing  the  product  of  quadra¬ 


ture  weight  and  the  function  value  Y(x^)  evaluated  at  the  collocation 


point  x^.  In  Equation  Dll  the  value  of  can  be  calculated  by  the 
following  equation: 


(2n  +  a  +  p  +  1)  C 


(a,P) 


W.  = 


^  x.  (1  -  x)  (x.)  y 

1  n  X 


(D12) 


where 


^  P^  (x)  (1  -  x)“  x^  dx 


(D13) 


Radau  quadrature  formula 

7.  In  Equation  Dll  the  end  points  Y(0)  and  Y(l)  are  not  required  to 
be  known.  If  either  one  of  the  end  points  is  given,  the  quadrature 
integration  can  be  obtained  more  accurately  by  using  the  Radau  formula. 

In  this  case  the  interior  quadrature  points  should  be  chosen  as 
if  Y(x=l)  is  given,  or  if  Y(x=0)  is  given.  The  Radau  formula 


D3 


can  integrate  a  polynomial  of  degree  up  to  2n  exactly  by  the  following 


equation; 


(1  -  x)“  Y(x)  =  Y(l)  +  1  W.  Y. 

i=l 


=  W  Y(0)  +  I  W.  Y. 
n+1  '  '  .  ,  11 

1=1 


The  Radau  quadrature  weights  in  Equation  D14  are  calculated  by  the 
following  equations; 

Including  *0~^‘ 


(DU) 


(2n  +  Of  +  p  +  2) 

-i  (■‘ill  (-i) 


1  ,  i#n+l 
l/(cf+l),  i=n+l 


(D15) 


Including  Xq=0,  but  not 


(2n  +  a  +  p  +  1) 

■  (I  -  X.)  (P<JJ  (X,)  )2 


1/(P+1),  i=0 

1,  i^O 


(D16) 


Computer  Programs  D1  and  D2 

8.  This  computer  program  is  for  computing  the  roots  of  the  Jacobi 

polynomial,  the  derivatives  of  each  polynomial  evaluated  at  the  colloca¬ 
tion  points,  and  the  quadrature  weights.  The  program  consists  of  two 
subroutines;  Subroutines  JCOBI  computes  the  zeros  of  ^nd  the 

first  derivatives  of  the  Jacobi  polynomial;  Subroutine  RADAU  determines 
the  integration  weights  at  collocation  points. 

9.  The  roots  of  the  Jacobi  polynomial  and  Radau  quadrature  weights 

for  selected  collocation  points  are  listed  in  Table  Dl.  An  example  of 

3  2 

calculating  the  integration  q/  x  dx,  using  the  roots  and  weights  found 
from  computer  program  Dl,  is  illustrated  in  computer  program  D2. 


''mV**' 
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COMPUTER  PROGRAM  D1 


this  program  is  used  to  compute  the  roots  and 

QUADRATURE  WEIGHTS  OF  JACOBI  POLYNOMIALS 
mmmm  DIFl,  DIF2,  AND  DIF3  ARE  THE  1ST.  2ND.  AND  3RD 

DERIVATIVES  OF  JACOBI  POLYNOMIALS  AT  COLLOCATION 
POINTS,  RESPECTIVELY 

ROOT(I)  ARE  ROOTS  OF  JACOBI  POLYNOMIALS 
WT(I)  ARE  QUADRATURE  WEIGHTS 
MMxw  N  IS  NUMBER  OF  COLOCATION  POINTS 

ALPHA  AND  BETA  ARE  THE  TWO  PARAMETERS  OF 
JACOBI  POLYNOMIAL 

IMPLICIT  REAL*8(A-H.  O-Z) 

DIMENSION  DIFl (11) ,DIF2( 1 1 ) ,DIF3( 1 1 ) ,ROOT( 1 1 ) ,WT( 1 1 ) 

1  FORMAT (215) 

3  FORMAT ( 1 H 1 . '  GEOMETRY :  0=PLANAR , 1 =CYLINDER , 2=SPHERE ' , 

&  //.*  NUMBER  OF  COLLOCATION  POINTS  =',15./,'  GEOMETRY 

&  FACTOR  =',15.///,  '  COLLOCATION  POINTS  IN  X*'*2  :  '  ./) 
6  F0RMAT(/.3(D16.8) ) 

200  READ(5. 1 .END=100)  N.IS 

N  IS  NUMBER  OF  COLLOCATION  POINTS 
FOR  PLANE  SHEET  IS=0,  FOR  CYLINDER  IS=1,  FOR  SPHERE  IS=2 

S  =  IS 

WRITE  (6,3)  N,  IS 
ALPHA  =  1.0 
BETA  =  (S-1)/2 

CALL  JCOBI  ( 1 1 ,N,0, 1 .ALPHA, BETA, DIFl ,DIF2,DIF3, ROOT) 
NT  =  N  +  1 

WRITE  (6,6)  (ROOT(I),I  =  1 ,NT) 

. FIND  QUADRATURE  WEIGHTS 

CALL  RADAU  ( 1 1 , N , 0 , 1 , 1 , 0 . DO , BETA , ROOT , DIF  1 , WT ) 

. TO  COMPUTE  TRUE  QUADRATURE  WEIGHTS  FOR  SLAB  GEOMETRY 

. THE  WEIGHTS  WT(I)  OBTAINED  FROM  SUBROUTINE  RADAU  NEED 

. TO  MULTIPLY  BY  "2.0- 

DO  90  I  =  1  ,  NT 
90  WT(I)  =  WT(I)'«2.0 
WRITE(6.55) 

55  FORMAT(/,3X, '  THE  QUADRATURE  WEIGHTS  WT(I)',/) 
WRITE(6,6)  (WT(I),  1=1 , NT) 

C 

CO  TO  200 
100  STOP 
END 


nnn  onnn  nnnn  nnnn  nnnnnnnnnnn 


COMPUTER  PROGRAM  D1  (CONTINUED) 


SUBROUTINE  JCOBI (ND ,N .NO , N 1 .AL.BE.FA.FB.FC.ROOT) 
SUBROUTINE  JCOBI  IS  USED  TO  COMPUTE  THE  ROOTS  OF 
JACOBI  POLYNOMIALS 

ND  IS  DIMENSION  OF  VECTORS,  N  IS  THE  DEGREE  OF 
JACOBI  POLYNOMIAL.  NO  DECIDES  WHETHER  X=0  IS  INCLUDES 
AS  AN  INTERPOLATION  POINT.  NO  MUST  BE  SET  EQUAL  TO 
1  (INCLUDING  X=0)  OR  0  (EXCLUDING  THIS  POINT) 

N1  IS  THE  SAME  AS  FOR  NO,  BUT  FOR  THE  POINT  X=1 
FA.FB.FC  ARE  1ST,  2ND,  AND  3RD  DERIVATIVES  OF 
JACOBI  POLYNOMIAL  AT  THE  NODES 

AL=ALPHA,  BE=BETA,  ROOT=ZEROS  OF  JACOBI  POLYNOMIAL 

IMPLICIT  REALi«8(A-H.0-Z) 

DIMENSION  FA(ND) ,FB(ND) .FC(ND) .ROOT(ND) 

THE  FIRST  STEP  IS  TO  CALCULATE  GN  &  HN 
HERE  FA(I)=GN,  FB(I)=HN 

AB=AL+BE 
AD=BE-AL 
AP=BE*AL 


. INITIAL  VALUE 

. FA( 1 )  IS  G1 ,  FB( 1 )  IS  H( 1 ) 

FA( I )=(AD/(AB+2. )+1 . )/2. 

FB( 1 )=0. 

IF  (N  .LT.  2)  GO  TO  15 

. THEN  CONFUTE  FA(2),  FB(2) ,  FA(3) ,  FB(3) .  ETC. 

. FA(I)=GN,  FB(I)=HN 

DO  10  I  =  2,N 
Z1  =I  -  1 
Z  =  AB  +  2.*Z1 

-  FA(I)  IS  GN.  FB(I)  IS  HN 

FA(I)=(AB*AD/Z/(Z+2. )+1 . )/2. 

IF(I  .NE.  2)  GO  TO  11 
FB(I)=(AB  +  AP  +  Z1 )/Z/Z/(Z+1 . ) 

GO  TO  10 
1 1  Z=Z*Z 

Y=Z1*(  AB  +  ZD 
Y=Y*(AP+Y) 

FB(I)=Y/Z/(Z-1 . ) 

10  CONTINUE 


COMPUTER  PROGRAM  D1  (CONTINUED) 


THE  SECOND  STEP 

ROOT  DETERMINATION  BY  NEWTON  METRHOD  WITH  SUPPRESION 
OF  PREVIOUSLY  DETERMINED  ROOTS 


X=0. 

DO  20  1=1 .N 
XD=0. 

XN=1  . 

XD1=0. 

XN1=0. 

DO  30  J=1 .N 

FA(J)=GN,  FB(J)=HN 

XP  IS  JACOBI  POLYNOMIAL,  AND  XP1  IS  THE  FIRST 
DERIVATIVE  OF  JACOBI  POLYNOMIAL 

XP=(FA(J)-X)*XN-FB( J)*XD 

XP1=(FA(J)-X)*XN1-FB(J)*XD1-XN 

XD=XN 

XD1=XN1 

XN  AND  XN1  ACCUMULATE  THE  FOREGOING  RESULTS 

XN=XP 
XN1=XP1 
ZC=1  . 

2=XN/XN 1 

IF(I  .EQ.  1  )  GO  TO  21 
DO  22  J  =  2,  I 

EXCLUDE  PREVIOUS  DETERMINED  ROOT:  ROOT(J-I) 
ZC=ZC-Z/ ( X-ROOT ( J- 1 ) ) 

Z=Z/ZC 
FIND  NEW  X 

x=x-z 

IF  (  DABS(Z)  .GT.  1 .D-09  )  GO  TO  25 
ROOT(I)=X 

NEW  STARTING  POINT  FOR  NEXT  ROOT 

X=X  +  .0001 
CONTINUE 

ADD  EVENTUAL  INTERPOLATION  POINTS  AT  X=0  OR  X=1 


n  n  n  n 


COMPUTER  PROGRAM  D1  (CONTINUED) 


NT=N+NO+N1 

IF (NO  .EO.  0)  GO  TO  35 
DO  31  1=1 .N 
J  =  N  +  1  -  I 
31  R00T(J+1)=R00T(J) 

R00T( 1 )=0. 

35  IF(N1  ,EQ.  1)  ROOT(NT)  =  1. 

. NOW  EVALUATE  DERIVATIVES  OF  POLYNOMIAL 

_  FA,  FB,  FC  ARE  1ST,  2ND,  AND  3RD  DERIVATIVES 

DO  40  1=1 , NT 
X=ROOT(I) 

FA(I)=1 . 

FB(I)=0. 

FC(I)=0. 

DO  40  J=1 ,NT 

IF  (J  .EQ.  I)  GO  TO  40 

Y=X  -  ROOT(J) 

FC(I)  =  FC(I)  *  Y  +  3.  «  FB(I) 

FB(I)  =  Y  *  FB(I)  +  2.  *  FA (I) 

FAC  I)  =  Y  *  FA (I) 

40  CONTINUE 
RETURN 
END 
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COMPUTER  PROGRAM  D1  (CONTINUED) 


SUBROUTINE  RADAU (ND .N .NO .N 1 , ID , AL .BE , ROD ,FA , V) 

C 

C . SUBROUTINE  RADAU  IS  USED  TO  EVALUATE  RADAU  OR  LOBATTO 

C  QUADRATURE  WEIGHTS  OF  JACOBI  POLYNOMIAL 

C . FAC  I)  COME  FROM  JCOBI,  V(I)  ARE  QUADRATURE  WEIGHTS 

C . ND,  N,  AND  NO  ARE  THE  SAME  AS  FOR  SUBROUTINE  JCOBI 

C . ID=1  GIVES  RADAU  QUADRATURE  WITH  X= 1 

C . ID=2  GIVES  RADAU  QUADRATURE  WITH  X=0 

C . ID=3  GIVES  LOBATTO  QUADRATURE  WITH  BOTH  POINTS 

C 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  RODCND) .FA(ND) ,V(ND) 

S=0. 

NT=N+N0+N1 
DO  40  1=1 .NT 
X=ROD(I) 

IF  (ID-2)  10,20,30 
C 

C  ID  =  1  GIVES  RADAU  QUADRATURE  WITH  X  =  1 

C  ID  =  2  GIVES  RADAU  QUDRATURE  WITH  X  =  0 

C  ID=3  GIVES  LOBATTO  QUADRATURE  WITH  BOTH  ENDPOINTS 

10  AX=X 

IF  (NO)  11,11,40 

11  AX=1./AX 
GO  TO  40 

20  AX=1.-X 

IF  (ND  21,21,40 

21  AX=1./AX 
GO  TO  40 

30  AX=1. 

40  V(I)=AX/FA(I)»<»<2 
IF  (ID-2)  41 ,42,41 

41  V(NT)=V(NT)/'(  1  .+AL) 

42  IF  (ID-2)  43,44,44 
44  V( 1 )=V( 1 )/( t .+BE) 

43  DO  50  1=1 , NT 
C 

C . S=V(  1)  +  V(2)  +  . 

C 

50  S=S+V(I) 

DO  60  1=1 .NT 
C 

C  V(I)  ARE  NORMALIZED  SO  THAT  THE  SUM  OF  V(I)  EQUALS  TO  ONE 
C***  NOTE  THE  QUADRATURE  WEIGHTS  COMPUTED  FROM  RADAU  ARE  NOT 
C  TRUE  WEIGHTS 
C 

60  V(I)=V(I)/S 

RETURN 
END 
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COMPUTER  PROGRAM  D2 


CALCyLATING  THE  INTEGRATION  OF  X  BY  USING  RADAU 
OUADI^TURE 

INTEGRATION  INTERVAL:  FROM  0  TO  3 


IMPLICIT  REAL»8(A-H,0-Z) 

DIMENSION  ROOTC 1 1 ) ,WT( 1 1 ) ,XX( 1 1 ) 

DATA  ROOT/0. 00533697,  0.04735263,  0.12782106, 

t  0 . 2399 1861,  0 . 374 1 3950 ,  0.51910186, 

^  0.66251286,  0.79221088,  0.89719596, 

Si  0.96855602,  1.00000000/ 

DATA  WT/0 . 29 1 69804 ,  0.28548098,  0.27317938, 

t  0.25505539,  0.23149529,  0.20300115, 

L  0 . 1 70 1 80 1 2 ,  0 . 1 3373 121,  0 . 09442893 , 

L  0.05309150,  0.00865801/ 


SUM  =  O.DO 
X  =  3. DO 


DO  10  I  = 
XXd)  =  X 
CONTINUE 


DSQRT(ROOT(I) ) 


DO  20  I  = 
SUM  =  SUM 
CONTINUE 


,  11 
WT(I) 


XXd)  *  XXd) 


SUM  =  3./2.«  SUM 


WRITE(6,1)  SUM 
FORMAT (/ 1 5X,  ’  rX  DX  = 


F10.6) 


SOLUTION;  j  X  DX  =  9.000000 


7**2 

w 


m 


Table  D1.  The  Roots  and  Quadrature  Weights  of 

Jacobi  Polynomial  for  10-points  Collocation 


Geometry :  Slab 


Collocation  Points  in 


ES; 

0 . 00533697 

0 . 04735263 

0.12782106 

Ei- 

0.23991861 

0.37413950 

0.51910186 

0.66251286 

0.79221088 

0.89719596 

• 

'r'i- 

0.96855602 

1 .00000000 

Si' 

Quadrature  Weights; 

0.291 69804 

0.28548098 

0.27317938 

0.25505539 

0.23149529 

0.203001 15 

0.17018012 

0. 13373121 

0.09442893 

0.05309150 

0.00865801 

vMvivivl. 


vv; 


APPENDIX  E:  INTEGRATION  OF  SEDIMENT  VOLUME 


1.  In  this  study,  Simpson’s  rule  is  used  for  the  calculation  of 
the  sedimenjt  volume  deposited.  Detailed  formulation  is  discussed  by 
Scheid  (1968).  Computation  procedure  and  computer  programs  are  briefly 
described  in  the  following. 

2.  The  volume  integration  of  a  two-variable  sediment  thickness 
function  f(x,y)  can  be  expressed  as 

fRbo  rb(x)  ^bg 

I  I  f(x,y)dydx  =  j  G(x)dx  (El) 

t  *0  0 

where  Rb^  is  the  dimensional  sediment  length  in  the  x-direction,  expressed 
as  the  Xj- direction;  and  b(x)  is  the  width  of  sediment  in  the  y-direction, 
expressed  as  the  x^-direction. 

By  setting  y^  =  0 

and  y  -  y  =  h,(x)  = 

J  J  n 


where  h^  is  the  discrete  distance  in  y-direction  and  n  is  the  number  of 
discretizations , 


^b(x) 

G(x)  =  /  f(x,y)dy 

0 


(  f(x,yQ)  +  4f(x,yj)  +  2f(x,y2)  + 


(E2) 


Again,  setting 


Xo  =  0 


and 


Rbf 

X.  -  X.  =  h  = 

1  1-1  2  m 


El 


JtkVl.lL 


where  h^  is  the  discrete  distance  in  x-direction  and  m  is  the  number  of 
discretizations , 


j  “  G(x)dx  =  [G(Xq)  + 


4G(Xj)  +  IGU^)  +  .... 

4G(x  )  +  G(x  )] 
m~  1  m 


Therefore 


f(x,y)dydx  =  ^ 


hi(Xo) 

“1 —  ^  + 


2f(Xo,y2) 


+  f(x. 


(E3) 


4h^(xp 


(f(Xi,yo)  +  4f(Xj,yj)  +  ZfCx^.y^)  +  -  +  f(Xj,y^)] 


2hj^(x2) 


3 -  (f(x2.yQ)  ^f(x2,yj)  +  2f(x2,y2)  + 


+  f(x2,yQ)] 


U(x^,yo)  +  +  2f(x_^,y2)  +  ....+  f(x^,yj] 


(E4) 


The  integration  of  the  foregoing  equation  is  illustrated  in 
The  numbers  in  Figure  El  represent  the  coefficients  of  f(x^. 
Equation  E4. 


Figure 
y^)  in 


El. 


Computer  Program  El 


3.  This  program  is  for  plotting  a  two-dimensional  sediment  thick¬ 
ness  for  a  given  set  of  parameters.  Using  the  equations  derived  in  PART 
V,  the  SAS/GRAPH  is  executed  for  displaying  the  plot.  The  longitudinal 
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Figure  El.  Integration  coefficients  of  Simpson's  rule 
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COMPUTER  PROGRAM  El 


DATA  ONE; 

UO= 1.0;  W0=0 . 05E-3 ; F=0 .001;  80=500 . ;  C0=300 . ;  E=0 . 075 ; 
H0=2.0;  AI1=0.948;  AI2=1.397;  AI3=0.76; 

TIME=2. »365. *86400. ;  DS=0.4; 

ROOT  1 =0 . 00533697 ;  R0OT2=0 . 04735263 ;  ROOT3=0 . 1 2782 1 06 
ROOT4=0 . 2399 1 86 1 ;  R00T5=0 . 374 1 3950 ;  ROOT6=0 . 5 1 91 0 1 86 
ROOT7=0. 6625 1286;  R0OT8=0. 79221 088;  ROOT9=0 . 897 1 9596 
ROOT  1 0=0 . 96855602 ;  ROOT 11  =  1. 00000000 ; 

WT 1=0. 29 169804;  UT2=0 . 28548098 ;  WT3=0 . 273 1 7938 ; 

WT4=0 . 25505539 ;  WT5=0 . 23 1 49529 ;  WT6=0 . 20300 115; 

WT7=0 .17018012;  WT8=0 . 1 3373 121;  WT9=0 . 09442893 ; 

WT 1 0=0 . 05309 150;  WT 1 1 =0 . 0086580 1 ; 

A=0.0001 ; 

SI  =  0.0;  S2=0.2;  S3=0.4;  S4=0.5;  S5=0.6;  S6=0.8; 

UCR  =  UO; 

W  =  BO*WO/UO/HO; 

AA  =  1.794*E; 

FF  =  F*BO/HO; 

D  =  F/A; 

DO  R  =  0.0  TO  20.0  BY  1.0; 

H  =  1.  +  A«BO/HO»R; 

CF  =  H0/A/'B0/(2.-D)«(  H**(2.-D)  -  1.)  +  0.5/AA; 

U  =  1 ./SQRT(2. »AA)  /  H»»D  /  SQRT(CF) ; 

B=  2. *AA«H** (D-1 )  ■  CF; 

HUl  =  1.  +  A*BO/HO«R*SORT(ROOT1 ); 

HU2  =  1.  +  A*BO/HO»R*SQRT(ROOT2) ; 

HU3  =  1.  +  A«BO/HO*R*SQRT(ROOT3) ; 

HU4  =  1,  +  A«BO/HO»R*SQRT(ROOT4) ; 

HU5  =  1.  +  A*BO/HO«R*SQRT(ROOT5) ; 

HU6  =  1.  +  A*BO/HO*R*SQRT(ROOT6) ; 

HU7  =  1.  +  A*BO/HO*R*SQRT(ROOT7) ; 

HU8  =  1.  +  A*BO/T10*R*SQRT(ROOT8)  ; 

HU9  =  1.  +  A«BO/HO*R*SQRT(ROOT9) ; 

HUIO  =  1.  +  A*BO/HO*R*SQRT(ROOT10) ; 

HU11  =  1.  +  A*BO/HO«R*SQRT(ROOT1 1 ) ; 

CFU1  =  H0/A/B0/(2.-D)*(  HU1*«(2.-D)  -  1.)  +  0.5/AA; 

CFU2  =  H0/A/B0/(2.-D)*(  HU2»»(2.-D)  -  1.)  +  0.5/AA; 

CFU3  =  H0/A/B0/(2.-D)*(  HU3»*(2.-D)  -  1.)  +  0.5/AA; 

CFU4  =  H0/A/B0/<2.-D)»(  HU4»»(2.-D)  -  1.)  +  0.5/AA; 

CFU5  =  H0/A/B0/(2.-D)*(  HU5*»(2.-D)  -  1.)  +  0.5/AA; 

CFU6  =  H0/A/B0/(2.-D)»(  HU6«»(2.-D)  -  1.)  +  0.5/AA; 

CFU7  =  H0/A/B0/(2.-D)«(  HU7«*(2.-D)  -  1.)  +  0.5/AA; 

CFU8  =  H0/A/B0/(2.-D)»(  HU8»*(2.-D)  -  1.)  +  0.5/AA; 

CFU9  =  H0/A/B0/(2.-D)*(  HU9»*(2.-D)  -  1.)  +  0.5/AA; 

CFU10  =  H0/A/B0/(2.-D)*(  HU10»«(2.-D)  -  1.)  +  0.5/AA; 

CFUII  =  H0/A/B0/(2.-D)»(  HUI1«»(2.-D)  -  1.)  +  0.5/AA; 

UU1  =  1./S0RT(  2.*AA)  /  HU  1 *»D/SQRT(CFU 1) ; 

UU2=1./SQRT(  2.*AA)  /  HU2*»D/SQRT(CFU2) ; 

UU3=1./SQRT(  2.*AA)  /  HU3*»D/SQRT(CFU3) ; 
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COMPUTER  PROGRAM  El  (CONTINUED) 


UU4=1./SQRT(  2.*AA)  /  HU4»«D/SQRT(CFU4) ; 

UU5=l,/SORT(  2.*AA)  /  HU5»»D/SQRT(CFU5) ; 

UU6=1>SQRT(  2.«AA)  /  HU6* » D/SORT ( CFU6 )  ; 

UU7=1./SQRT(  2.«AA)  /  HU7««D/SQRT(CFU7) ; 

UU8=1./SQRT(  2.»AA)  /  HUS* • D/SORT ( CFU8 ) ; 

UU9=1./SQRT(  2.*AA)  /  HU9*«D/SQRT(CFU9) ; 

UU10=1  ./SQRT(  2.*AA)  /  HU  1 0 »"•  D/SORT (CFU  1  0)  ; 

UU1 1  =  1 ./SQRT(  2.»AA)  /  HU  11 **D/SQRT(CFU 1 1) ; 

G1=  -AI2*W/AI1»'R/2.  ; 

G2=  AI3*W/AI1/UCR/UCR*R/2. ; 

SUMl  =  WT1/HU1/UU1  +  WT2/HU2/UU2  +  WT3/HU3/UU3 
+  WT4/HU4/UU4  +  WT5/HU5/UU5  +  WT6/HU6/UU6 
+  WT7/HU7/UU7  +  WT8/HU8/UU8  +  WT9/HU9/UU9 
+  UT 1  O/HU 1  0/UU 1 0  +  UT1  1/HU1  1/UUn  ; 

SUM2  =  WT1«UU1/HU1  +  WT2«UU2/'HU2  +  WT3«UU3/HU3 
+  WT4«UU4/HU4  +  WT5*UU5/HU5  +  WT6»UU6/HU6 
+  WT7*UU7/HU7  +  WT8*UU8/HU8  +  WT9*UU9/HU9 
+  WT10*UU 10/HU  10  +  WT1 1*UU1 1/HUl 1 ; 

C  =  EXP(  G1»SUM1  +  G2*SUM2  -  LOG(H*U«B)  ); 

CA1  =  CO«SQRT( 1 .-S1»S1 )/EXP(Sl«S1/2. )*C; 

CA2  =  COi'SQRTC  1  .-S2*S2)/EXP<S2»S2/2.  )»C; 

CA3  =  CO*SQRT(  1  .-S3'«S3)/EXP(S3»S3/2.  )«C; 

CA4  =  CO-SQRTC 1 .-S4*S4)/EXP(S4*S4/2. )*C; 

CAS  =  CO*SQRT( 1 .-S5*S5)/EXP(S5*S5/2. )«C; 

CA6  =  CO*SQRT( 1 .-S6*S6)/EXP(S6«S6/2. )*C; 

CCA1=<1.  -  S1»S1)/EXP(S1*S1 ); 

CCA2=(1.  -  S2»S2)/EXP(S2»S2); 

CCA3=(1.  -  S3*S3)/EXP<S3*S3) ; 

CCA4=(1.  -  S4*S4)/EXP(S4»S4); 

CCA5=(1.  -  S5*S5)/EXPCS5*S5); 

CCA6=<1.  -  S6*S6)/EXP(S6»S6) ; 

UA1=U0*CCA1»U; 

UA2=U0*CCA2«U; 

UA3=U0»CCA3»U; 

UA4=U0»'CCA4*U; 

UA5=U0»CCA5*U; 

UA6=UC«'CCA6*U; 

CDl  =  W0*100.»CA1«1 ,0E-6«<1 .  -  UA 1 «UA 1/UCR/UCR) »TIME 

CD2  =  WO* 100. *CA2« 1 .0E-6*( 1 .  -  UA2*UA2/UCR/UCR) *TIME 

CD3  =  W0*I00.*CA3*1 .0E-6*( 1 .  -  UA3*UA3/UCR/UCR) *TIME 

CD4  =  W0*100.*CA4*l.0E-6*(1.  -  UA4-UA4/UCR/UCR) *TIME 

CDS  =  WO* 100. *CAS« 1 .0E-6»( 1 .  -  UAS*UAS/UCR/UCR ) *TIME 

CD6  =  W0*100.*CA6*1 .0E-6*<1 .  -  UA6«UA6/UCR/UCR ) *TIME 

TH1  =  CD1/DS; 

TH2  =  CD2/DS; 

TH3  =  CD3/DS; 

TH4  =  CD4/DS; 

TH5  =  CDS/DS; 

TH6  =  CD6/DS; 
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COMPUTER  PROGRAM  El  (CONTINUED) 


OUTPUT; END; 

TITLE  1  .F=TRIPLEX  .H=2  PLOT  OF  SEDIMENT  THICKNESS; 

TITLE2  .R=TRIPLEX  .H=l  :IN  CM,  TWO  YEARS  OF  DEPOSITION; 
TITLES  .F=TRIPLEX  .H=1  CASE  4:  B0=500.,  U0=1.0,  W0=0.05; 
TITLE4  .F=TRIPLEX  .H=1  H0=2.0,  F=0.001; 

TITLES  .F=TRIPLEX  .H=1  S  =  0.0:  DIAMOND,  S  =  0.2:  STAR; 

TITLES  .F=TRIPLEX  .H=l  S  =  0.4:  PLUS,  S  =  0.5:  SQUARE 

TITLE7  .F=TRIPLEX  .H=l  S  =  0.6:  TRIANGLE.  S  =  0.8:  HASH; 

PROC  GPLOT; 

PLOT  TH1»R  TH2*R  TH3»R  TH4*R  TH5*R  TH6-R  /  OVERLAY; 

SYMBOL  1  I=SPLINE  V=DIAMOND; 

SYMB0L2  I=SPLINE  V=STAR; 

SYMBOLS  I=SPLINE  V=PLUS; 

SYMB0L4  I=SPLINE  V=SQUARE; 

SYMBOLS  I=SPLINE  V=TRI ANGLE; 

SYMBOLS  I=SPLINE  V=HASH; 


distance  RV  is  taken  corresponding  to  the  tnaximum  sediment  extent  from 
the  plot  (for  the  sake  of  computation,  RV  is  used  as  corresponding  to 
0.5  cm  sediment  thickness). 


4.  For  Case  4,  linearly  varying  depth  with  entrainment,  the 
following  hominal  values  are  given: 


bg  =  500  m 
Uq  =  1.0  m/sec 
f  =  0.001 
Cq  =  300  ppm 
=  400  kg/m^ 


Wq  =  0.05  mm/ sec 

hg  =  2.0m 
a  =  0.0001 
e  =  0.075 


The  roots  and  quadrature  weights  of  collocation  for  integration,  taken 
from  Appendix  D  as  inputs,  are  shown  in  Table  El. 

5.  The  resulting  plot  of  sediment  thickness  of  this  example  is 
shown  in  Figure  E2;  the  longitudinal  distance  RV  is  19.7.  The  plots  for 
other  bifurcation  levels  can  be  generated  in  a  like  manner  by  supplying 
the  appropriate  values  for  b^. 


Computer  Program  E2 


6.  Using  the  longitudinal  distance  RV  from  the  plot  of  Computer 
Program  El  and  other  given  parameters  as  inputs,  this  FORTRAN  computer 
program  evaluates  total  volume  of  sediments  VOL  by  Simpson's  rule  for  the 
integration.  The  variable  names  in  the  computer  program,  written  as 
similarly  as  possible  to  the  expressions  shown  in  PART  V,  are  not  ex¬ 
plained  in  detail  here. 

7.  For  Case  4,  the  RV's  and  bg's  are  given  as  following: 

RV  =  19.7,  51.9,  130.0,  229.0,  and  387.5 
bg  =  500.,  166.7,  55.6,  27.8,  and  13.9 

The  calculated  volume  of  sediments  are  shown  in  Table  E2. 
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COMPUTER  PROGRAM  E2 


C 

C  EVALUATE  TOTAL  SEDIMENT  FOR  GIVEN  BO 

C  CASE  4,  REAL  DIMENSION 

C 

IMPLICIT  REAL»8(A-H.  0-2) 

C 

DATA  DS/0.4D0/.  UO/1 .  OODO/'. HO/2.  ODO/.CO/SOO .  DO/ 

C 

TIME  =2.0 

100  READ(5.51 .END=999)  RV,  BO 

51  F0RMAT(2F10.4) 

C 

CALL  INTEG(RV,BO.TIME,VOL) 

C 

c 

WRITE(6,1)  RV.BO.VOL 

1  FORMAT( IHl ,///'  RV  =',F10.4,/.'  BO  =',F9.4,'  (M) ' ./, 

&  '  VOLUME  OF  SEDIMENT  VOL  (M«»3)  ='.D16.6) 

C 

GO  TO  100 
999  STOP 
END 

SUBROUTINE  INTEG(RO,BO,TIME, VOL) 

IMPLICIT  REAL»8(A-H.  0-Z) 

DIMENSION  HI ( 1 1 ) ,THICK( 1 1 , 1 1 ) .A( 11 ) .B( 11 ) ,C( 1 1 ) 

C 

C  Y  IS  X2  COORDINATE,  X  IS  XI  COORDINATE 

C 

RBO  =  RO-BO 
C 

E  =  0.075 
F  =  0.001 
HO  =  2.0 
AA  =  1 .794»E 
FF  =  F-BO/HO 
SA  =  0. ID-3 
D  =  F/SA 
M  =  5 
N  =  5 
C 

DO  20  J  =  1 .  11 

XI  =  (J  -  n/lO.  *  RBO 

H2  =  RBO  /M/2. 

C 

C  Bl=S*B(>:i),  NOTE  :  S=1.0 

C  FIND  JET  WIDTH  AT  VARIOUS  XI 

C 

R  =  Xl/BO 

HH  =  1 .  +  SA«BO/HO«R 
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COMPUTER  PROGRAM  E2  (CONTINUED) 


CF  =  H0/SA/B0/(2.-D)  *  (HH*» (2. -D)-1 . )  +  0.5/AA 
B1  =  BO  *  2.  »  AA  «  HH»»(D-1.)  «  CF 
HI (J)  =  B1  /N  /2. 

'*• 

DO  1 0  I  =  1 ,  11 

X2  =  (I  ~  1)/10.  *  B1 

CONFUTE  SEDIMENT  THICKNESS  AT  (XI ,X2) 

CALL  BED(B0,TIME.X1 ,X2.B1 ,THK) 

THICKd.J)  =  THK 
1 0  CONTINUE 
20  CONTINUE 


THUS  ALL  VALUES  OF  THICKd.J)  HAVE  BEEN  OBTAINED 


START  CALCULATING  SEDIMENT  VOLUME 


FIND  END  POINT  ON  THE  R  COORDINATE  FOR  SIMPSON'S 
RULE 

DO  30  I  =  1 ,  11 
30  Ad)  =  THICKd,  1  ) 

H  =  Hid) 

CALL  SIMPS(N.H,A,SUM) 

SUMTI  =  H2/3.  «  SUM 

DO  40  I  =  1 ,  11 
40  Bd)  =  THICKd,  1  1  ) 

H  =  HI ( 1 1 ) 

CALL  SIMPS(N,H.B.SUM) 

SUMT2  =  H2/3.  «  SUM 


COMPUTE  ALL  VALUES  BETWEEN  TWO  END  POINTS 

SUMT3  =0.0 
IFLAG  =  I 
DO  50  J  =  2,  10 
H  =  HI (J) 

DO  60  I  =  1 ,  11 
Cd)  =  THICKd.J) 

60  CONTINUE 

CALL  SIMPS(N,H,C.SUM) 

IF (IFLAG  .EO.  1 )  GO  TO  111 
IF (IFLAG  .EQ.  2)  GO  TO  222 
1 1 1  SUM  =  H2/3.  ■4.  »  SUM 
SUMT3  =  SUMT3  +  SUM 
IFLAG  =  2 
GO  TO  50 
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COMPUTER  PROGRAM  E2  ( CONTINUED > 

222  SUM  =  H2/3.  «  2.  ■  SUM 
SUMT3  =  SUMT3  +  SUM 
IFLAG  =  1 
50  CONTINUE 

SUMTT  =  SUMT1  +  SUMT2  +  SUMT3 

:  VOL  IS  VOLUME  OF  SEDIMENT  IN  (M»«3) 

VOL  =  SUMTT 

RETURN 

END 


SUBROUTINE  SIMPSCN.HI .FUNC.SUMM) 

IMPLICIT  REAL»8<A-H.  0-Z) 

DIMENSION  FUNCdn 

USE  SIMPSON'S  RULE  FOR  INTEGRATION 
SUMEND  IS  SUM  OF  ALL  FUNC(I)  FOR  EVEN  I 
SUMMID  IS  SUM  OF  ALL  FUNC(I)  FOR  ODD  I 
HI  IS  STEP  SIZE 
FUNC  IS  INTEGRAND 


INITIALIZE  PARAMETERS 

SUMEND  =0,0 
SUMMID  =0.0 

EVALUATE  SUMEND  AND  SUMMID 

DO  1  K  =  1  ,  N 
KI  =  2mK  -  1 
K2  =  2*iK 

SUMEND  =  SUMEND  +  FUNC(Kl) 

SUMMID  =  SUMMID  +  FUNC(K2) 

IT  =  2«N  +  1 

SUMM  =  (2.0*(SUMEND-FUNC( I ) )  +  4.0 "SUMMID  +  FUNCC 1 ) 
&  +  FUNCCIT))  *  Hl/S.O 
RETURN 


SUBROUTINE  BED(B0 .TIME.XI ,X2,BX1 .THICK) 

IMPLICIT  REAL*8(A-H.  0-Z) 

DIMENSION  ROOT( 1 1 ) ,WT< 11 ) 

DATA  HO/2 .  ODO/ ,  A/0 .0001  DO/ ,  AI 1  /0 . 948D0/ ,  AI2/ 1 . 397D0/ 
&  AI3/0.76D0/,  DS/0.4D0/,  E/0.075D0/,  UO/l.ODO/, 

&  F/O.OOIDO/,  W0/0.05D-3/,C0/300.D0/,  UCR/l.ODO/ 

DATA  ROOT  /0.005337D0,  0.047353D0,  0.1 2782 IDO, 


COMPUTER  PROGRAM  E2  (CONTINUED) 


&  0.239919D0,  0.374139D0,  0.519102D0, 

B.  0.662513D0,  0.79221100,  0.897196D0, 

i  0.968556D0.  l.OOOOOODO/ 

DATA  UT  /  0.291698D0. 
i  0.255056D0. 

i  0.170180D0. 

i  0.05309  IDO, 

CUCR  =  UCR/UO 
U  =  BO«WO/UO/HO 
D  =  F/A 
AA  =  1 .794«E 
FF  =  F»BO/HO 

R  =  Xl/BO 
S  =  X2/BX1 
HH  =1.  +  A«BO/'HO*R 

G1  =  -AI2»W/AIl»R/2. 

G2  =  AI3*U/'AIl/CUCR/CUCR»R/2. 

AG1  =  HO/A/BO/ ( 2. -D)  »  (HH»* (2.-D)-1 . )  +  0.5/AA 

U  =  1 ./DSQRTC2.*AA)/  HH«"D  /DSQRT(AGl) 

B  =  2. »AA«HH«*(D-1 )  «  AGl 

SUM1  =0.0 
SUM2  =0.0 

DO  1 0  I  =  1 .  11 

HU  =  1.  +  A»BO/HO»R»DSQRT(ROOT(I) ) 

CFU  =  HO/A/BO/ ( 2. -D)  «  (HU»'"(2.-D)-1  .  )  +  O.S/AA 
UU  =  1 ./DSQRT(2."AA)/  HU*«D  /DSQRT(CFU) 

SUM1  =  SUM1  +  WTCD/HU/UU 
SUM2  =  SUM2  +  WT(I)/HU«UU 
CONTINUE 

C  =  DEXP(  Gl-SUMl  +  G2-SUM2  -  DLOG ( HH»U»B) ) 

Cl  =  CO«DSQRT( 1 .-S"S)/DEXP(S»S/2. )«C 
CC1  =  (1.  -  S»3)/DEXP(S«S) 

U1  =  U0»CC1»U 

THICK  IS  THICKNESS  OF  SEDIMENT  IN  METER 

THICK  =  W0»C1»1 .0D-6«( 1 .-U1-U1/UCR/UCR) 

,  »TIME«365*86400./DS 

RETURN 

END 


0.285481D0,  0.273180D0, 
0.231495D0,  0.203001D0, 
0 . 1 3373 IDO,  0 . 094429D0 , 
0.008658D0  / 


Computer  Program  E3 


8.  Using  50  percent  of  the  total  volume  of  sediments  calculated 
from  Computer  Program  E2  and  other  given  parameters  as  inputs,  this  com¬ 
puter  program  searches  for  the  longitudinal  distance  RA  of  the  deposited 
sediments  by  the  bisection  method. 

9.  For  Case  4,  after  two  years  of  deposition,  the  total  sediments 

6  3 

volume  VOL  is  59.58  x  10  m"^  for  the  1st  generation.  The  normalized 

longitudinal  distance  RA,  for  a  given  volume  of  total  sediments,  is  found 
to  be  9.36  after  11  iterations.  A  sample  output  of  this  case  is  shown  in 
Table  E3. 


Computer  Program  E4 

10.  This  computer  program  is  used  for  searching  the  time-step 
required  to  fill  the  known  volume  (which  is  50  percent  of  the  total  vol¬ 
ume)  of  sediments  to  an  average  thickness.  Using  the  longitudinal  dis¬ 
tance  RV  of  the  total  volume  (calculated  from  Computer  Program  E2)  and 
RA  of  the  50  percent  total  volume  (calculated  from  Computer  Program  E3) 
as  inputs,  the  time-step  is  obtained  by  the  method  of  bisection  method. 
Also  the  dimensional  distance  of  Xj^  (length),  x^  (width),  total  area  AREA, 
and  total  volume  VOL  of  the  deposited  sediments  are  calculated. 

11.  For  Case  4,  first  generation, 

RV  =  19.7 

RA  =  9.36 

A  listing  of  output  for  this  case  is  shown  in  Table  E4. 
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COMPUTER  PROGRAM  E3 

SEARCH  FOR  R  OF  DEPOSITION  CORRESPONDING  TO  50% 

OF  TOTAL  VOLUME 

CASE  4,  REAL  DIMENSION 

USING  BISECTION  METHOD  FOR  SEARCH 

IMPLICIT  REAL«8(A-H,  0-Z) 

DATA  DS/0.4D0/',  U0/1.0D0/,  H0/2.0D0/,  C0/300.D0/ 

XN  AND  XP  ARE  THE  SEARCH  RANGE  OF  R 
VOL1  IS  THE  TOTAL  VOLUME 

RV  IS  LONGITUDINAL  DISTANCE  OF  TOTAL  VOLUME 

100  READ(5,51 ,END=999)  XN,  XP,  BO,  RV,  VOL1 

51  F0RMAT<4F10.4,D20.4) 

WRITE(6.3n  BO.  RV.  VOLl 

31  FORMAT(//.'  BO  =  •.F10.4,*  (M) V, '  RV  =' .F10.4./, 

&  '  INITIAL  VOLUME  VOLl  =  ',D16.6.’  (M»»3) ' ) 

ITER  =  0 

12  XM  =  (XN  +  XPV2. 

ITER  =  ITER  +  1 

IF(  DABS(XP-XN)  .LT.  0.01  )  GO  TO  10 

DIF  IF  THE  DIFERENCE  BETWEEN  EXPECTED  AND 
COMPUTED  VOLUMES 

TIME  =2.0 
VOL2  =  VOLl  »  0.5 
R  =  XM 

CALL  INTEG(BO,R,TIME,VOL) 

DIF  =  VOL  -  V0L2 

WRITE(6,1)  ITER.  DIF,  VOL,  XM 
1  FORMATC/.'  NO.  OF  ITERATION  =',I5,/, 

&'  ERROR  OF  VOLUME  (M*»3)  ='.D16.6,/, 

&'  VOLUME  OF  SEDIMENT  (M-xS)  =',D16.6,/, 

&'  RA  =' ,F12.4,/) 

IF(DIF)  11,  10,  13 
11  XN  =  XM 
CO  TO  12 

13  XP  =  XM 
GO  TO  12 

10  WRITE(6.21)  XM 

21  FORMATC/,’  mmmmmm  RA  =',  F12.5,/) 
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COMPUTER  PROGRAM  E3  (CONTINUED) 

GO  TO  100 
999  STOP 
END 

SUBROUTINE  INTEG(B0.R0. TIME, VOL) 

IMPLICIT  REAL»8(A-H,  0-Z) 

DIMENSION  HI ( 1 1 ) .THICKC 11 , 1 1 ) , A( 1 1 ) ,B( 11 ) ,C( 1 1 ) 

Y  IS  X2  COORDINATE,  X  IS  XI  COORDINATE 

RBO  =  R0>«B0 

E  =  0.075 
F  =  0.001 
HO  =  2.0 
AA  =  1 .794«E 
FF  =  F«B0/H0 
SA  =  0. ID-3 
D  =  F/SA 
M  =  5 
N  =  5 

DO  20  J  =  1 .  11 
XI  =  ( J  -  1 )/10.  *  RBO 
H2  =  RBO  /M/2. 

B1=S*B(X1),  NOTE  :  S=1.0 
FIND  JET  WIDTH  AT  VARIOUS  XI 

R  =  X1/B0 

HH  =  1 .  +  SA*B0/H0«R 

CF  =  H0/SA/B0/(2.-D)*(HH»»(2.-D)-1 . )  +  0.5/AA 
B1  =  80*2. ■AA*HH«*(D-1 . )*CF 
HHJ)  =  B1  /N  /2. 

DO  1 0  I  =  1  ,  11 
X2=  (I  -  1)/10.  «<  B1 

CONFUTE  SEDIMENT  THICKNESS  AT  (XI. X2) 

CALL  BED(B0,TIME,X1 .X2,B1 ,THK) 

THICKd.J)  =  THK 
10  CONTINUE 
20  CONTINUE 

THUS  ALL  VALUES  OF  THICKd.J)  HAVE  BEEN  OBTAINED 
START  CALCULATING  TOTAL  VOLUME 

FIND  END  POINT  ON  THE  R  COORDINATE  FOR  SIMPSON'S 
RULE 


nnnnnn  n  n  n  n  nnn 


COMPUTER  PROGRAM  E3  (CONTINUED) 


DO  30  I  =  1 ,  11 
30  Ad)  =  THICKd.l) 

H  =  HKD 

CALL  SIMPS<N,H.A.SUM) 

SUMT1  =  H2/3.  «  SUM 

DO  40  I  =  1 ,  11 
40  B(I)  =  THICKCI, 1 1 ) 

H  =  HI ( 1 1 ) 

CALL  SIMPS(N.H.B,SUM) 

SUMT2  =  H2/'3.  ■  SUM 

COMPUTE  ALL  VALUES  BETWEEN  TWO  END  POINTS 

SUMT3  =0.0 
IFLAG  =  1 
DO  50  J  =  2,  10 
H  =  HI <J) 

DO  60  I  =  1 ,  11 
C(I)  =  THICKCI, J) 

60  CONTINUE 

CALL  SIMPS(N,H.C,SUM) 

IF( IFLAG  .EQ.  1 )  GO  TO  1 1 1 
IF (IFLAG  .EQ.  2)  GO  TO  222 
111  SUM  =  H2/3.  »4.  »  SUM 
SUMT3  =  SUMT3  +  SUM 
IFLAG  =  2 
GO  TO  50 

222  SUM  =  H2/3.  «  2.  «  SUM 
SUMT3  =  SUMT3  +  SUM 
IFLAG  =  1 
50  CONTINUE 

SUMTT  =  SUMT1  +  SUMT2  +  SUMT3 

VOL  IS  VOLUME  OF  SEDIMENT  IN  (M«»3) 

VOL  =  SUMTT 

RETURN 

END 

SUBROUTINE  SIMPSCN ,H1 .FUNC.SUMM) 

IMPLICIT  REAL*8(A-H.  0-Z) 

DIMENSION  FUNC( 1 1 ) 

USE  SIMPSON’S  RULE  FOR  INTEGRATION 
SUMEND  IS  SUM  OF  ALL  FUNC(I)  FOR  EVEN  I 
SUMMID  IS  SUM  OF  ALL  FUNC(I)  FOR  ODD  I 
HI  IS  STEP  SIZE 
FUNC  IS  INTEGRAND 
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COMPUTER  PROGRAM  E3  (CONTINUED) 


INITIALIZE  PARAMETERS 

SUMEND  =0.0 
SUMMID  =  0.0 

EVALUATE  SUMEND  AND  SUMMID 

DO  1  K  =  1  ,  N 
K1  =  2»K  -  1 
K2  =  2«K 

SUMEND  =  SUMEND  +  FUNC(KI) 

1  SUMMID  =  SUMMID  +  FUNC(K2) 

IT  =  2*N  +  1 

SUMM  =  (2.0*(SUMEND-FUNC( I ) )  +  4.0*SUMMID  +  FUNC( 1 ) 

&  +  FUNC(IT) )  «  HI/3.0 
RETURN 
END 

C========================================================= 

SUBROUTINE  BED ( BO .TIME , X 1 ,X2.BX1 .THICK) 

IMPLICIT  REAL»«8(A-H.  O-Z) 

DIMENSION  ROOT( 1 1 ) .WT( 1 1 ) 

C 

DATA  H0/2.0DO/.A/0.0001DO/,AI 1/0.948D0/.AI2/1 .397D0/ 
&  AI3/0.76D0/.  DS/0.4D0/,  E/0.075D0/.  UO/I.ODO/. 

&  F/0.001D0/.  WO/0.05D-3/.CO/300.DO/.  UCR/ 1 . ODO/ 

C 

DATA  ROOT  /0.005337DO.  0.047353DO.  0.1 2782 IDO. 

&  0.239919D0.  0.374139D0.  0.5191 02D0. 

&  0.662513D0,  0.79221  IDO.  0.897196D0. 

&  0.968556DO.  1.000000D0  / 

C 

DATA  WT  /  0.291698D0. 

&  0.255056D0. 

&  0.170180DO. 

&  0.05309  IDO. 

CUCR  =  UCR/UO 
W  =  BO-WO/UO/HO 
D  =  F/A 
AA  =  1 .794-E 
FF  =  F-BO/HO 
C 

R  =  XI /BO 
S  =  X2/BX1 
HH  =  1 .  +  AoBO/HOoR 
C 

G1=  -AI2-W/AI 1 »R/2. 

G2  =  AI3»W/AI  1./UCR/UCR«R/2. 

C 

AG1  =  H0/A/B0/(2.-D)*(HH«»(2.-D)-1 . )  +  0.5/AA 


0.28548  IDO,  0.273180D0. 
0.231495D0,  0.20300  IDO, 
0 .  1 3373  IDO,  0 . 094429D0 , 
0.008658D0  / 
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COMPUTER  PROGRAM  E3  (CONCLUDED) 

C 

U=  1 ./DSQRT(2. *AA)/HH**D/DSQR;(AG1 ) 

C 

B  =  2 . *AA*HH* * ( D- 1) *AG 1 
C 

SUM1  =0.0 
SUM2  =0.0 
C 

DO  1 0  I  =  1 ,  11 

HU  =  1.  +  A*BO/HO*R*DSQRT(ROOT(I) ) 

CPU  =  H0/A/B0/(2.-D)*(HU**(2.-D)-1 . )  +  0 . 5/AA 
UU  =  1 ./DSQRT^2. *AA)/HU**D/DSQRT(CFU) 

SUM1  =  SUM1  +  WT(I)/HU/UU 
SUM2  =  SUMP  ^  WKD/HU^UU 
10  CONTINUE 
C 

C  =  DEXP(G1*SUM1  +  G2*SUM2  -  DLOG ( HH»U»B ) ) 

Cl  =  C0*DSQRT( 1 .-S*S)/DEXP(S*S/2. )*C 
CC1  =  (1.  -  S*S)/DEXP(S*S) 

U1  =  U0»CC1*U 

THICK  IS  THICKNESS  OF  SEDIMENT  IN  METERS 

THICK  =  W0*C1*1 .0D-6*( 1 .-U1*U1/UCR/UCR) 

&  *TIME*365. *86400. /DS 

RETURN 
END 

C  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =:^^..  =  =  =  =  =  =  =  =:  =  ===  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  : 
//GO. SYS IN  DD  * 

5.0  10.0  500.0  19.7  0 . 5958D08 

// 


COMPUTER  PROGRAM  E4 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

100 

51 

c 


31 


C 


12 


C 

C 

C 

C 


c 

c 

c 


c 


1 


1 1 
13 
10 


SEARCH  FOR  TIME  OF  DEPOSITION  FOR  GIVEN  RV  (OF  VOLUME) 
AND  RA  (OF  AREA).  AND  AVERAGE  THICKNESS  OF  DEPOSITION 
CASE  4.  REAL  DIMENSION 

USING  BI-SECTIONAL  METHOD  FOR  SEARCH 

IMPLICIT  REAL*8(A-H,  0-Z) 

XN  AND  XP  ARE  THE  SEARCH  RANGE  OF  SEDIMENTATION  TIME 
RV  IS  R  FOR  VOLUME.  RA  IS  R  FOR  AREA 

READ(5.51 .END=999)RV.RA.XN.XP,B0 
F0RMAT(5F10.2) 

URITE(6.31)  BO.  RV.  RA 

FORMAT( IHl .///. '  BO  =  '.F10.4,'  (M) ’ 

8.  •  RV  =  •.F10.4./.*  RA  =  '  ,F10.4) 


ITER  =  0 

XM  =  (XN  +  XP)/2. 

ITER  =  ITER  +  1 

IF(  DABS(XP-XN)  .LT.  O^OIDO  )  GO  TO  10 

DIF  IF  THE  DIFERENCE  BETWEEN  EXPECTED  AND 
COMPUTED  AVERAGE  SEDIMENT  THICKNESS 

TIME  «  XM 

CALL  INTEG(BO.RV.TIME.VOL) 

UNIT  OF  AREA  =  M««2 

CALL  INTEG2(B0.RA.AREA.X2) 

A  =  0. ID-3 
HO  =  2.0D0 

DIF  =  VOL/AREA  -  (HO  +  0.5  *  A  »  RA  «  BO) 

WRITE(6.1)  ITER,  DIF,  VOL,  AREA,  XM 
FORMAT(/,'  NO.  OF  ITERATION  =',15./, 

&'  ERROR  OF  BED  THICKNESS  (M)  =' ,D16.6,/, 

&'  VOLUME  OF  SEDIMENT  (M»»3)  =' ,D16.6,/, 

8,'  AREA  OF  SEDIMENT  (M««2)  =' ,D16.6,/, 

&'  TIME  (YEARS)  =',F12.4) 

IF(DIF)  11,  10,  13 

XN  =  XM 

GO  TO  12 

XP  =  XM 

GO  TO  12 

XI  =  RA«BO 

WRITE(6,21)  XI.  X2,  BO.  XM 
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COMPUTER  PROGRAM  E4  (CONTINUED) 

FORMATC/,*  XI  =  ’.FIO-A,/,'  X2  =  ',F10.4,/, 
&  '  BO  =  ',F10.4,/, 

&  '  time  (YEARS)  =*,  F12.5) 

GO  TO  100 

STOP 

END 


SUBROUTINE  INTEG (BO ,R0 .TIME, VOL) 

IMPLICIT  REAL»<8(A-H.  O-Z) 

DIMENSION  HI ( 1 1 ) ,THICK( 11 , 1 1 ) . A( 11 ) ,B( 1 1 ) .C( 11 ) 

Y  IS  X2  COORDINATE,  X  IS  XI  COORDINATE 

RBO  =  R0»B0 

E  =  0.075 
F  =  0.001 
HO  =  2.0 
AA  =  1 .794«E 
FF  =  F*B0/H0 
SA  =  0. ID-3 
D=F/SA 
M  =  5 
N  =  5 

DO  20  J  =  1 ,  11 

XI  =  (J  -  1 )/10.  «  RBO 

H2  =  RBO  /M/2. 

B1=S*B(X1),  NOTE  :  S=1.0 
FIND  JET  WIDTH  AT  VARIOUS  XI 

R  =  Xl/BO 

HH  =  1 .  +  SA»BO/HO»R 

CF  =  H0/SA/B0/(2.-D)»(HH>«*(2.-D)-l  .  )+0.5/AA 
B1=B0»2.  ■'AA*HH»»  CD- 1  .  )*CF 
H  1  ( J )  =  B 1  /N  /2 . 

DO  1 0  I  =  1 ,  11 

X2  =  (I  -  1 )/10.  ■  B1 

CONFUTE  SEDIMENT  THICKNESS  AT  (XI ,X2) 

CALL  BED ( BO . TIME . X 1 , X2 , B 1 , THK ) 

THICK(I.J)  =  THK 

CONTINUE 

CONTINUE 

THUS  ALL  VALUES  OF  THICK ( I, J)  HAVE  BEEN  OBTAINED 
START  CALCULATING  TOTAL  VOLUME 


COMPUTER  PROGRAM  E4  (CONTINUED) 

FIND  END  POINT  ON  THE  R  COORDINATE  FOR  SIMPSON'S 
RULE 

DO  30  I  =  1 .  11 
30  A(I)  =  THICKd.l) 

H  =  Hl(l) 

CALL  SIMPS (N.H. A. SUM) 

SUMTl  =  H2/3.  ■  SUM 

DO  40  I  =  1  ,  11 
40  B(I)  =  THICKd.  1  1  ) 

H  =  HI ( 1  1  ) 

CALL  SIMPS(N.H.B,SUM) 

SUMT2  =  H2/3.  ■  SUM 

COMPUTE  ALL  VALUES  BETWEEN  TWO  END  POINTS 

SUMT3  =0.0 
IFLAG  =  t 
DO  50  J  =  2.  10 
H  =  HI <J) 

DO  60  I  =  1 ,  11 
Cd)  =  THICKd, J) 

60  CONTINUE 

CALL  SIMPS(N,H,C,SUM) 

IFdFLAG  .EQ.  1  )  GO  TO  1 1  1 
IF (IFLAG  .EQ.  2)  GO  TO  222 
111  SUM  =  H2/3.  »4.  «  SUM 
SUMT3  =  SUMT3  +  SUM 
IFLAG  =  2 
GO  TO  50 

222  SUM  =  H2/3.  ■  2.  *  SUM 
SUMT3  =  SUMT3  +  SUM 
IFLAG  =  1 
50  CONTINUE 

SUMTT  =  SUMTl  +  SUMT2  +  SUMT3 

VOL  IS  VOLUME  OF  SEDIMENT  IN  (M»«3) 

VOL  =  SUMTT 
RETURN 


SUBROUTINE  INTEG2(B0 
IMPLICIT  REAL«8(A-H, 

.R0.AREA.X2) 

0-Z) 

A  =  0. lD-3 

F  =  0.00  IDO 

Y'  -■ 


COMPUTER  PROGR_AM  E4  (CONTINUED) 

D  =  F/A 
S  =  1 .0 
E  =  0.075 
F  =  0.001 
HO  =  2.0 
AA  =  1 .794*E 
FF  =  F»BO/HO 


UNIT  OF  AREA  =  M«*2 

CK1  =  H0»»2/(2. *(A»B0)»*2*(2.-D) ) 

CK2  =  H0/’(D«A»B0)  »  <0.5/AA  -  H0/(A«B0)/(2. -D) ) 

HH  =  1 .  +  A  «  BO/HO  »  RO 

CF  =  HO/A/BO/ ( 2. -D)  *  (HH  »»  (2.-D)  -  1.)  +  O.B/AA 
X2  =  BO  ■  2.«  AA  »  HH*»(D-1.)  *  CF 
AREA  =  (2.  »<AA«B0*«2)  «  (CK1»(1.  +  A»BO/HO»RO )  «*2 
&  +  CK2  »  (1.  +  A*BO/'HO»RO)*«D  -  CK1  -  CK2) 

RETURN 
END 

SUBROUTINE  SIMPS ( N , H 1 , FUNC , SUMM ) 

IMPLICIT  REAL*8(A-H,  0-2) 

DIMENSION  FUNC(II) 

USE  SIMPSON'S  RULE  FOR  INTEGRATION 
SUMEND  IS  SUM  OF  ALL  FUNC(I)  FOR  EVEN  I 
SUMMID  IS  SUM  OF  ALL  FUNC(I)  FOR  ODD  I 
HI  IS  STEP  SIZE 
FUNC  IS  INTEGRAND 

INITIALIZE  PARAMETERS 

SUMEND  =0.0 
SUMMID  =0.0 

EVALUATE  SUMEND  AND  SUMMID 

DO  1  K  =  1  ,  N 
K1  =  2*K  -  1 
K2  =  2XK 

SUMEND  =  SUMEND  +  FUNC(KI) 

1  SUMMID  =  SUMMID  +  FUNC(K2) 

IT  =  2»N  +  1 

SUMM  =  (2.0*(SUMEND-FUNC( 1 ) )  +  4.0*SUMMID  +  FUNC( 1 ) 
&  +  FUNC(IT))  «  HI/3.0 
RETURN 
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COMPUTER  PROGRAM  E4  (CONTINUED) 

SUBROUTINE  BED (BO. TIME. XI ,X2.BX1 .THICK) 

IMPLICIT  REAL-8 (A-H.  O-Z) 

DIMENSION  ROOT( 1 1 ) ,WT( 11 ) 

C 

DATA  HO/2 . ODO/ . A/0 . 000 1 DO/ , AI 1 /O . 948D0/ , AI2/ 1 . 397D0/ . 
&  AI3/0.76D0/,  DS/0.4D0/,  E/0.075D0/,  UO/l.ODO/, 

&  F/O.OOIDO/,  W0/0.05D-3/,C0/300.D0/,  UCR/l.ODO/ 

C 

DATA  ROOT  /0.005337D0.  0.047353D0.  0.1 2782 IDO. 

0.239919D0,  0.374139D0,  0.519102D0, 

0 . 6625 1 3D0 ,  0 . 7922 1 1  DO ,  0 . 897 1 96D0 , 
0.968556D0,  1.000000D0  / 


C 


c 


c 

c 

c 

c 


c 


10 


c 


WT  /  0.291698D0. 
0.255056D0. 
0. 170180D0. 
0.05309 IDO, 
CUCR  =  UCR/UO 
U  =  BO-UO/UO/HO 
D  =  F/A 
AA  =  1 .794-E 
FF  =  F-BO/HO 


0.28548 IDO,  0.273180D0, 
0.231495D0,  0.203001D0, 
0 . 1 3373 IDO,  0 . 094429D0 , 
0.008658D0  / 


R  =  X1/B0 

S  =  X2/BX1 

HH  =  1 .  +  A-BO/HO-R 

Gl=  -AI2-W/AI1-R/2. 

G2=  AI3-W/AI1/CUCR/CUCR-R/2. 


AGI  =  H0/A/B0/(2.-D)*(HH»*(2.-D)-l . )  +  0.5/AA 


U  =  1 ./DSQRT(2.»AA)/HH**D/DSQRT(AG1 ) 


B=  2. *AA*HH-»(D-1 )-AG1 


SUM1  =0,0 
SUM2  =0.0 


DO  1 0  I  =  1 ,  11 

HU  =  1.  +  A*<BO/HO»R*DSORT(ROOT(I) ) 

CFU  =  H0/A/B0/(2.-D)»(HU»*(2.-D)-l . )  +  0.5/AA 
UU  =  1 ./DSQRT(2.*AA)/HU««D/DS0RT(CFU) 

SUMl  =  SUMl  +  WT(I)/HU/UU 
SUM2  =  SUM2  +  WT(I)/HU*UU 
CONTINUE 

C  =  DEXP(  Gl-SUMl  +  G2-SUM2  -  DLOG(HH*U»B) ) 

Cl  =  CO-DSQRT(  1  .-S'*S)/DEXP(S'«S/2.  )*C 
CC1  =  (1.  -  S*S)/DEXP(S»S) 

U1  =  UO-CCl-U 


E24 


c 

c 
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COMPUTER  PROGRAM  E4  (CONCLUDED) 

THICK  IS  THICKNESS  OF  SEDIMENT  IN  METERS 


THICK 

RETURN 

END 


& 


W0*C1 *  1 .0D-6*( 1 .-U1/UCR/UCR) 
*TIME*365. *86400. /DS 


//GO.SYSIN  DD  * 

19.7  9.36 

// 
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